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ABSTRACT 


A  numerical  code,  USPOTF2,  has  been  formulated  to  solve  for  the  potential  flow 
for  two  airfoils  executing  unsteady  motions  in  an  inviscid  incompressible  flow  medium. 
This  code  is  an  extension  of  an  existing  code  U2DIIF,  which  does  the  same  calculations 
for  the  single  airfoil  case.  The  technique  uses  the  well  known  Panel  Methods  for  steady 
flow  and  extends  it  to  unsteady  flow  by  introducing  a  wake  model  which  creates  a  non¬ 
linear  problem  due  to  the  continuous  shedding  of  vortices  into  the  trailing  wake.  The 
presence  of  the  second  airfoil  introduces  a  set  of  non-linear  coupled  equations  for  the 
Kutta  condition.  Numerous  case-runs  are  presented  to  illustrate  the  capability  of  the 
code.  The  case  of  the  step  change  in  angle  of  attack  is  compared  with  Giesing's  work. 
All  other  case-runs  are  illustrated  together  with  the  results  for  the  single  airfoil  case. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research  may  not 
have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made,  within 
the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic  er¬ 
rors,  they  cannot  be  considered  validated.  Any  application  of  these  programs  without 
additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 


A.  GENERAL 

In  this  paper  ,  a  numerical  method  is  formulated  to  solve  for  the  flow  about  two 
two-dimensional  airfoils  which  are  arbitrary-  located  and  are  performing  an  arbitrary  time 
dependent  motion  in  an  inviscid  incompressible  fluid.  The  original  work  by  Teng  [Ref. 
1]  is  for  a  single  arbitrarily  defined  airfoil  in  the  same  potential  flow  condition.  The  ex¬ 
tension  of  Teng's  code  to  two  bodies  is  considered  here.  Where  possible,  the  same  con¬ 
vention  and  notation  are  adopted.  As  for  the  single  airfoil  case,  all  velocities  are 
non-dimensionalised  with  respect  to  the  free  stream  and  all  lengths  with  respect  to  the 
chord  length. 

A  new  subroutine  (SUBROUTINE  NEWPOS)  is  added  to  transform  either  of  the 
tw  o  local  coordinate  systems  to  the  global  coordinate  system.  A  modification  is  entered 
into  the  treatment  of  the  Kutta  condition  to  make  it  consistent  with  the  unsteady  flow 
Kutta  condition  treatment.  Subroutine  COFISH  is  deleted  and  its  role  is  taken  over  by 
Subroutine  COEF  which  normally  performs  the  formulation  for  the  flow  tangency  con¬ 
dition  for  the  unsteady  case.  A  more  accurate  method  is  also  introduced  to  obtain  the 
velocity  potential.  The  reader  is  referred  to  the  w’ork  of  Krainer  (Ref.  2]  for  further 
improvement  of  the  original  code. 

This  documentation  is  set  up  as  for  the  original  documentation  in  order  that  it  will 
be  easier  to  follow-  and  cross-referenced.  While  it  is  the  intent  of  the  author  to  keep  this 
thesis  as  complete  as  possible,  the  reader  is  wreil  advised  to  review  Reference  1  for  the 
w'ork  involved  in  the  single  airfoil  case;  no  special  effort  will  be  made  to  reproduce  it 
here. 

B.  BASIC  THEORY  AND  APPROACH 
1.  Steady  Flow  Problem 

The  treatment  for  the  two  airfoils  case  in  steady  incompressible  flow  follow's 
closely  to  the  single  airfoil  case.  The  governing  equation,  as  for  the  single  airfoil  case  , 
follows  from  the  Conservation  of  Mass  and  the  Condition  of  Irrotational  Flow. 

The  continuity  equation  of  an  incompressible  fluid  ( div  V  -  0)  and  the  condi¬ 
tion  of  irrotational  flow'  {curl  V  -  0)  leads  to  the  well  known  Laplace  Equation. 

div(grad<£)  *  0  (1.1) 
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with  <t>  denoting  the  disturbance  potential  for  the  velocity.  This  is  seen  to  be  the  classic 
Neumann  problem  of  potential  theory'  with  the  usual  problem  of  defining  the  boundary 
conditions.  The  boundary  conditions  for  the  disturbance  potential  4>  are  that  its  gradi¬ 
ent  normal  to  the  surface  be  equal  to  the  normal  velocity  of  the  surface  and  that  its 
gradient  vanish  at  infinity;  that  is 

V<t>  •«  —  Vr*n  on  S  (1.2) 

lim  V^(  P)  -*  0  (1.3) 

P—ao 

where  V,  is  the  resultant  velocity  of  a  point  on  the  body  as  seen  from  the  inertia  frame 
of  reference,  n  is  the  outward  unit  vector  normal  to  the  body,  S  is  the  body  surface  and 
P  represents  a  general  point.  Equation  1.2  holds  for  both  airfoils  i.e.  on  both  surfaces. 

The  pressure  coefficient  is  obtained  through  the  Bernoulli's  equation  which  de¬ 
rives  from  the  Momentum  equation. 

The  approach  adopted  is  associated  with  Hess  and  Smith  [Ref.  3]  who  devised 
the  popularly  known  PANEL  method  in  the  early  sixties.  In  words,  the  boundary  or 
airfoil  surfaces  S,  and  S,  about  which  the  flow  is  to  be  computed  is  approximated  by  a 
large  number  of  surface  elements  whose  characteristic  dimensions  are  small  compared 
to  those  of  the  body.  Over  each  surface  element,  a  uniform  source  distribution  and  a 
uniform  vorticity  distribution  is  placed.  The  source  strength  ($,)  varies  from  element  to 
element,  while  the  vortex  strength  (y,)  is  the  same  for  all  elements  in  the  same  airfoil  but 
is  different  across  the  airfoil.  The  singularity  strengths  are  determined  from  the  flow 
tangency  condition  on  both  body  surfaces  and  the  two  Kutta  conditions  at  both  trailing 
edges.  With  the  determination  of  the  singularity  strengths,  the  relevant  aerodynamic 
data  can  then  be  subsequently  computed. 

2.  Unsteady  Flow  Problem 

The  unsteady  problem  is  similar  to  the  steady  flow  problem  in  that  they  both 
have  the  same  governing  equation  viz.  the  Laplace  equation,  and  that  for  both  problems, 
the  pressure  and  velocity  are  decoupled  so  that  the  velocity  and  pressure  calculation  can 
be  computed  separately  and  consecutively. 

This  problem  differs  from  the  steady  flow  in  that  another  model  is  required  to 
simulate  the  continuous  shedding  of  vorticity  into  the  trailing  wake.  The  existence  of  a 
vortex  sheet  behind  the  airfoil  can  be  explained  by  the  Helmholtz  theorem  w'hich  is 
basically  a  statement  of  the  Conservation  of  Vorticity.  This  requires  that  any  change  in 
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the  circulation  around  the  airfoil  must  be  matched  by  an  equal  and  opposite  vortex 
somewhere  in  the  fiowfield.  The  presence  of  the  countervortices  provides  the  flow  with 
a  kind  of  a  memory  in  that  the  flow  at  a  particular  time  is  affected  by  the  bound  circu¬ 
lation  of  the  past.  It  is  this  non-linearity  that  distinguishes  the  numerical  technique  from 
the  simple  steady  flow  problem  of  solving  N  linear  equations  in  N  unknowns. 

The  solution  technique  requires  an  iterative  type  solution.  The  present  ap¬ 
proach  follows  closely  the  original  panel  method  of  Hess  and  Smith  as  described  in  the 
steady  flow  development,  while  with  regard  to  the  modelling  of  the  wake,  it  adopts  the 
procedure  advocated  by  Basu  and  Hancock  [Ref.  4).  A  uniform  source  distribution  {q)k 
and  a  uniform  vorticity  distribution  [y(/)l*  as  for  steady  flow  is  placed  on  each  panel  at 
time  t„  where  j  denotes  the  panel  number  and  /  the  airfoil  number.  The  wake  consists 
of  a  single  vorticity  panel  attached  as  an  additional  element  on  each  airfoil  through 
which  the  vorticities  are  shed  into  the  respective  wake  and  a  series  of  point  vortices 
which  are  being  converted  downstream  with  the  fluid.  A  uniform  vorticity  distribution 
of  strength  ('/.(/))*  is  placed  on  the  wake  panel  of  each  airfoil.  This  panel  is  further 
characterised  by  its  length  A(/)*  and  its  inclination  ( 0^,,, )» with  respect  to  the  respective 
local  frame  of  reference.  After  each  time  step,  the  vorticity  of  the  wake  panel  is  con¬ 
centrated  into  a  single  point  vortex  and  converted  downstream.  Simultaneously  a  new 
wake  panel  is  formed.  The  downstream  wake  of  point  vortices  is  thus  formed  by  the 
shed  vorticity  of  previous  time  steps. 

C.  SCOPE 

Chapter  II  extends  the  original  code  to  handle  the  steady  flow  problem  for  two 
airfoils  set  at  different  relative  distances  and  angles  of  attack. 

Chapter  III  deals  with  the  unsteady  problem  for  the  two  airfoils  system.  It  intro¬ 
duces  a  new  subroutine  and  a  more  accurate  method  of  calculating  the  velocity  poten¬ 
tial.  The  Kutta  condition  for  the  twro  airfoils  system  is  specially  treated  as  its  unique 
problem  of  a  non-linear  coupled  system  is  not  seen  in  the  single  airfoil  case. 

Chapter  IV  describes  the  computer  program,  its  essential  capabilities  and  limita¬ 
tions,  its  associated  subroutines,  its  input  requirements  and  its  associated  output  print¬ 
out. 

Chapter  V  presents  the  results  of  some  case-runs.  Of  interest  is  a  comparison  case 
of  a  step  change  in  angle  of  attack  (AOA)  with  Giesing  [Ref.  5]  for  the  same  airfoil 
undergoing  an  impulsive  start  at  the  same  AOA.  Case-runs  will  also  be  run  with  both 
airfoils  at  large  distances  apart  to  compare  with  the  single  airfoil  case.  In  addition,  ex- 
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ample  cases  for  which  no  comparisons  exist  are  given,  to  exhibit  the  capability  of  the 
method. 

Finally,  Chapter  VI  concludes  with  future  development  efforts  and  the  application 
potential  of  this  numerical  method. 
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II.  STEADY  FLOW  PROBLEM  FORMULATION  FOR  TWO  AIRFOILS 

A.  GENERAL 

The  modification  work  for  the  steady  flow  is  straightforward.  The  revised  program 
allows  for  two  arbitrarily  defined  airfoils  placed  at  an  arbitrary  distance  set  at  different 
angles  of  attack.  For  reason  of  simplicity,  the  number  of  panels  and  nodes  and  the  pivot 
location  are  set  to  be  the  same  for  both  airfoils. 

B.  FRAMES  OF  REFERENCE 

Three  frames  of  reference  are  involved  in  the  two  two-dimensional  airfoils'  case  in 
steady  flow.  These  are  three  inertia  frames  of  references  as  indicated  in  Figure  1. 

The  first  inertia  frame  of  reference  (also  known  as  global  frame  of  reference)  is  set 
at  the  pivot  position  of  the  first  airfoil  with  the  X-axis  pointing  in  the  direction  of  the 
free-stream  velocity.  The  two  other  inertia  frames  of  references,  henceforth,  will  be 
known  as  two  frozen  local  frames  of  referencel  and  x2y2)  are  fixed  respectively  to 
each  airfoil  with  the  x-axis  coinciding  with  the  chord  line  originating  from  the  respective 
leading  edge.  The  two  local  frames  of  reference  are  set  apart  by  XShift  and  YShift  on 
the  global  frame  of  reference.  In  steady  flow',  the  fluid  velocity  and  pressure  depend  only- 
on  the  spatial  coordinates  (X,Y)  and  not  on  time. 

The  two  airfoils  are  defined  in  the  local  coordinate  system  as  input  data  for  sim¬ 
plicity  and  are  then  transformed  to  the  global  coordinate  system  through  a  knowledge 
of  the  relative  positions  of  the  2  airfoils'  pivots  positions  and  the  respective  local  angles 
of  attack. 

C.  STEADY  FLOW  PANEL  METHODS 

1.  Definition  of  nodes  and  panels 

Each  airfoil  surface  is  divided  into  n  straight  line  segments  called  panels  by 
(n+ 1)  arbitrarily  chosen  points  called  nodes.  The  numbering  sequence  begins  with 
panel  1  on  the  lower  surface  at  the  first  airfoil  trailing  edge  and  proceeds  clockwise 
around  the  airfoil  contour  so  that  the  last  panel  on  airfoil  1  ends  on  the  upper  surface 
at  the  trailing  edge.  This  numbering  sequence  then  proceeds  in  a  similar  fashion  for  the 
second  airfoil  (  See  Figure  2).  As  with  the  single  airfoil  case,  the  numbering  sequence 

1  For  the  steady  case,  the  notation  'frozen'  will  be  dropped 


dictates  that  the  airfoil  body  always  lies  on  the  right  hand  side  of  the  i*  panel  as  one 
proceeds  from  the  i*  node  to  the  (i  +  1)"*  Also  the  1"  and  (n  +  l)1*  nodes  coincide  at  the 
first  airfoil  trailing  edge  and  the  (n  +  2)*  and  (2 n  +  2)rt  coincide  at  the  second  airfoil 
trailing  edge.  This  numbering  system  facilitates  the  definitions  of  the  unit  normal  vector 
rt,  and  the  unit  tangential  vector  t,  for  all  panels  with  nt  being  directed  outward  from  the 
body  into  the  flow  and  i,  directed  from  the  i*  node  to  the  (/  +  l)rt  node.  The  numbering 
of  the  panel  system  is  somewhat  complicated  by  the  fact  that  a  continuous  panel  num¬ 
bering  sequence  across  the  airfoil  is  desired.  This  procedure  leads  to  the  peculiar  panel 
numbering  behaviour  seen  in  Figure  2  where  for  the  first  airfoil,  the  panel  lies  between 
f*  and  (/  + 1)“*  nodes  while  for  the  second  airfoil,  the  panel  lies  between 
(/+  1)“  and  (/  +  2)*  nodes2. 

2.  Distribution  of  Singularities 

Over  each  surface  element  of  the  two  airfoils,  a  uniform  source  distribition  and 
a  uniform  vorticity  distribution  is  placed.  The  source  strength  q,  varies  from  element  to 
element  for  each  airfoil  w  hile  the  vorticity  strength  y,  remains  the  same  for  all  elements 
in  the  same  airfoil  but  is  different  across  the  airfoil.  This  choice  of  singularities  follows 
closely  the  original  panel  method  of  Hess  and  Smith.  It  automatically  satisfies  the 
Laplace  Equation  (which  is  the  governing  equation  for  the  inviscid  incompressible  flow) 
and  the  boundary  condition  at  the  far  field  (oo).  In  addition,  as  the  Laplace  Equation 
is  a  linear  homogeneous  second  order  partial  differential  equation,  an  overall  compli¬ 
cated  flow  field  can  be  built  up  by  the  combination  of  simple  flows  with  the  condition 
that  the  appropriate  boundary  condition  on  the  airfoil  be  satisfied  accurately. 

For  our  case,  the  overall  (low  field  (represented  by  the  velocity  potential  <D)  can 
be  built  up  by  three  simple  flows.  Writing  this  in  terms  of  the  respective  local  frame  of 
reference, 

,y)  -  4>UX  >y)  +  W*  >y)  +  <£»(*  >y)  (2.1) 

where  4>m{x  ,y)  is  the  potential  of  the  onset  flow, 

<t>ao(x ,  y)  *  FoJjcCos  a(/)  +  j/Sin  a (/)]  (2.2) 


2  For  the  code  6,  has  been  defined  out  of  the  above  convention  with  $„  i*  1,2..  ji  as  per  panel 
numbers  for  the  first  airfoil,  but  with  reserved  for  the  wake  element  of  the  first  airfoil  and  8„ 
i*n  +  2n+  3...2n  +  1  for  the  second  airfoil  with  reserved  for  the  wake  element  of  the  second 
airfoil  for  unsteady  flow. 
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2n+2 

"■C 

n+2 


Panel  j 

Source  distribution  »  qj 
Vorticity  distribution  «  ^ 


Note:  1.  Nodal  points  defined  by  i  =  1,2  ... 

and  i  *  n+2  . . . 


n+1  for  first  airfoil 
2n+2  for  second  airfoil. 


2. Panel  number  defined  by  j  *  1,2  . 
_  _  and  j  «  n+1  . 


n  for  first  airfoil 
2n  for  second  airfoil. 


Figure  2.  Panel  Methods  Representation  for  Stead)  Flow. 


<t>,  is  the  velocity  potential  of  a  source  of  distribution  q(s)  per  unit  length. 


*,« 


In  rds 


<t>,  is  the  velocity  potential  of  a  vorticity  distribution  y  (s)  per  unit  length 


(2.3) 


At  this  point,  the  disturbance  potential,  <f>,  is  introduced,  which  is  defined  to  be  the  sum 
of  the  potential  due  to  the  source  and  vorticity  distribution. 


*-*,  +  *»  (2-5) 

Equation  (2.1)  can  then  be  read  as 

4>-*oo  +  *  (2-6) 

The  convenience  of  defining  the  above  allows  for  the  total  velocity  vector  to  be  viewed 
as  two  components  viz.  the  onset  and  the  induced  velocity  due  to  the  disturbance  po¬ 
tential.  The  total  velocity  is  thus  : 

^  total m 

=  +  V<t>  (2.7) 

-?oo  +  V<f> 


It  is  in  the  introduction  of  the  disturbance  potential  that  leads  us  to  the  concept  of  in¬ 
fluence  coefficient  which  will  be  elaborated  in  a  later  section. 

The  pressure  coefficient  can  be  obtained  from  Bernoulli's  Equation  which  is 
derived  from  the  Conservation  of  Momentum. 

v  -  (2.8) 

fpvl  l  K“  > 
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3.  Boundary  Condition 

As  in  the  single  airfoil  case,  the  boundary  conditions  to  be  satisfied  include  the 
flow  tangency  conditions  and  the  Kutta  condition.  The  flow  tangencv  conditions  are 
satisfied  at  the  exterior  mid-points  (control  points).  The  normal  velocity  is  taken  with 
respect  to  the  respective  local  frame  of  reference  (for  consistency  with  unsteady  flow 
notation  to  be  introduced  later): 

(K^-O,  i- 1,2 .../!,«+  1  ...2n  (2.9) 

where  (K"),  is  the  normal  component  of  the  total  velocity*. 

The  Kutta  condition  postulates  that  the  pressure  on  the  upper  and  lower  sur¬ 
face  at  the  trailing  edge  of  each  airfoil  be  equal.  For  steady  potential  flow,  equal  pres¬ 
sure  implies  equal  tangential  velocity  in  the  downstream  direction  at  the  first  and  last 
panel  of  each  airfoil  viz  the  Bernoulli  equation.  With  our  definition  of  the  tangential 
vector  we  then  have 

{V\  -~(V%  1st  airfoil 

(V')n^  -  -(Y'hn  2nd  airfoil  (2.10) 

As  with  the  single  airfoil,  equations  2.9,  and  2.10  lead  to  a  linear  system  of 
(2n+2)  simultaneous  equations.  With  K"and  V'  expressed  explicitly  in  terms  of  q,  (j  - 
1,  2, ...  n,  n+ 1...  2n)  and  y,  (1®  1, 2)  we  have  2n + 2  unknowns  in  (2n+  2)  system  of  linear 
simultaneous  equations  which  can  be  easily  solved. 

D.  INFLUENCE  COEFFICIENT 
1.  Concept  of  Influence  Coefficient 

As  introduced  earlier,  this  important  concept  of  influence  coeflicient  results 
from  the  presence  of  the  disturbance  potential  which  follows  from  the  presence  of  the 
singularities.  Formally,  an  influence  coefficient  is  defined  as  the  velocity  induced  at  a 
Held  point  by  a  unit  strength  singularity  (be  it  a  point  singularity  or  a  distributed 
singularity)  placed  anywhere  within  the  flow  field.  Recall  that  equations  2.9  and  2.10 
require  the  computation  of  the  normal  and  tangential  velocity  components  at  all  the  el¬ 
ements  control  points.  The  normal  components  of  velocities  are  essential  in  satisfying 

3  This  follows  Giesing's  notation  where  the  velocity  with  respect  to  the  airfoil  frame  is  denoted 
as  total  velocity.  While  there  is  no  misunderstanding  in  the  steady  flow,  the  notation  actually  refers 
to  the  relative  velocity  with  respect  to  the  airfoil  moving  frame  of  reference  for  the  unsteady  flow. 


the  flow  tangency  conditions,  while  the  tangential  components  of  velocities  are  necessary 
for  satisfying  the  Kutta  condition  as  well  as  computing  the  pressure  distribution. 

2.  Notation  for  Influence  Coefficient 

The  notation  used  in  this  documentation  will  be  as  for  the  notation  for  the 
single  airfoil  case  with  a  slight  modification  to  take  into  account  the  effects  of  the  the 
second  airfoil.  As  the  influence  coefficients  are  related  to  the  geometry  of  the  airfoil  and 
their  relative  positions,  it  must,  of  necessity  be  computed  with  respect  to  the  global 
frame  of  reference  for  the  two  airfoils'  case. 

For  steady  flow,  the  following  influence  coefficients  are  defined. 

•  A" :  normal  velocity  component  induced  at  the  i*  control  point  by  unit  strength 
source  distribution  on  the  j*  panel. 

•  A'v :  tangential  velocity  component  induced  at  the  /**  control  point  by  unit  strength 
source  distribution  on  the  j*  panel. 

•  B i  :  normal  velocity  component  induced  at  the  i*  control  point  by  unit  strength 
vorticity  distribution  on  the  j*  panel. 

•  B *  :  tangential  velocity  component  induced  at  the  i*  control  point  by  unit  strength 
vorticity  distribution  on  the  /*  panel. 

where  i  and  j  denotes  panel  numbers  and  are  defined  as  : 

i*m  I,  2, ...  n,  n+  1  ...  2n 

j*m  1,  2, ...  n,  n+  1  ...  2n 

3.  Computation  of  Influence  Coefficient 

A  single  source  located  on  the  j*  panel  in  the  local  frame  of  reference  (see  Figure 
3)  induces  a  total  velocity  of 


(2.11) 


where  the  component  perpendicular  to  the  inducing  panel  is 


V" 


I  h 
2nr  r 


1 


2 *  h2  +  s2 


(2.12) 


and  the  component  parallel  to  the  inducing  panel  is 


Vs 


1 


2nr 


•  (-  —  )  « - — 

v  r  ’  In 


h2  +  s2 


(2.13) 


For  a  distributed  source  panel  on  the  j*  panel  in  the  local  frame  of  reference,  we  have 


K  ~  — £*  FTA\ 


(2.14) 


where 


FT  AN  *  tan 


(xm,  -  -  j^,)  -  (ym^y^xm,  -  *,+,) 

(*"•/  -  XjXxm,  -  ac;+1)  +  (ym,  -  .ty)(ym,  -  >)+l) 


Jf».  -  l/2(jf/  + Jf<+1) 


ym,  -  1/2  (y,  +.yJ+,) 


(2.15) 


^  -  f5jf7*)<* - 4^  FLOG 

Jst 


(2-16) 


where 


FLOG  -  In 


(xm,~xJ+y)2  +  Qi, -.v^,)2 
(arm,  -  */  +  (ym,  -y/ 


(2-17) 


Transforming  the  above  to  global  frame  of  reference,  we  have 


Ah  -  l^Cos (0,-dj)  -  Ipinte,- 6j) 


(2.18) 


A[  «  lpin(0,  -  6j)  +  F*Cos(e,  -  Bj) 


(2.19) 


Bl - (lyCos^  -  Gj)  -  F’Sin(0,  -  0,)) 

By  -  ^Sin^-fy)  -  l^Cos (0,-0,) 


(2.20) 


(2.21) 


For  the  Steady  Flow  code,  we  define 


AAN(1,J)  -  ^  -  Bj 


(2.22) 


flfl.W)  -  -Wj  -  B” 


(2.23) 
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E.  NUMERICAL  SOLUTION  SCHEME 
1.  Rewriting  the  Boundary  Condition 

Though  it  is  not  critical  for  the  steady  flow  problem  to  adopt  any  particular 
frame  of  reference,  we  will  adopt  the  local  frame  of  reference  in  satisfying  the  boundary 
condition  for  consistency  with  the  treatment  used  in  the  unsteady  case. 

Even  though  the  influence  coefficients  are  defined  in  terms  of  the  global  frame 
of  reference,  the  computation  of  the  normal  and  tangential  velocities  of  each  panel  is 
independent  of  the  coordinate  system  used.  Thus  the  normal  and  tangential  velocities 
obtained  with  the  global  coordinates  system  would  be  the  same  as  those  obtained  for  the 
local  coordinate  system.  Using  the  local  coordinate  system,  the  flow  tangency  condition 
of  equation  2.9  is  defined  as  follows : 


n  2n  n  2n 

(A  -  + Y + y(i)YB‘j + y{2)Y  B"y + v°° SinHf) - d,]  “  °  (i24) 

Jm\  J-n+l  Jm\  Jmit+1 

where 

/  —  1,2  ...  n  ;/-  1, 

*  n  +  1 ...  2n  ;  **  2 

The  Kutta  condition  of  eqn  2.10,  in  terms  of  influence  coefficients,  becomes  for  airfoil 

1: 

it  2  n  n  2  it 

Ya^  +  Y  + y{])YBxj + v(2)  Y Bxj + v°°  c°sca(i)  ~  ^  * 

jm\  7*n+ 1  Jm  1  y«n+l 


n  2n  n  2/1 

-  +YJA'rfJ+ y^Y8* + y(2)  YB'/+y°°  c°s[a(i)  _ 

U  y— I  Jrnfl+l  Jm  1  Jmn+\  ' 


(2.25) 


and  for  airfoil  2: 


it  2  it  it  2n 

An+\jQ]  +  y(f)^^^/«+iy  +  y(2)  Y^j  +  V~  Uos[a(2) 


>1 


Jmlt+l 


J- 1 


Jmn+\ 
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n  In  n  In 

g^lnjR]  +  ^lnj* 1)  +  y(l)^^2n,/  +  y(2)  ^  ^2, 

y—n+i  y»i  y*«+i 


2/iy  +  ^  Cos[a(2)  -  e2„]  l  (2.26) 


2.  Solving  for  the  Strengths  of  Source  end  Vorticity  Distribution 

Equations  2.24,  2.25  and  2.26  can  be  written  as  a  set  of  2n  +  2  linear  simultane¬ 
ous  equations  with  2n  +  2  unknowns  ( qJtj *■  1,2  ...  In  and  ylt  /*  1,2)  and  solved  as  for 
the  single  airfoil  case.  However,  to  make  the  routine  consistent  with  the  unsteady  flow 
case,  the  following  method  is  adopted. 

The  flow  tangency  condition  can  be  rewritten  explicitly  for  the  qf  in  terms  of 
y(l),  y(2)  and  the  free  stream  constant  term;  that  is 


■ 

" 

- 

■ 

*1,1 

<*1.2 

*U 

• 

• 

• 

*1,2/i 

ft 

*l,2n+l 

*1,271+2 

*2.1 

*2,2 

*2.3 
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• 

*2,2" 

ft 

*2 ,2/1+1 

*2.2/t+2 

*2 

*3.1 

*3.2 

*3.3 
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• 

*3,2/i 

ft 

*3,2/1+ 1 

*3 .2/1+2 

*3 
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y(2)  + 
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• 

• 

• 

• 

• 

• 

• 

• 

• 

*2/i,l 

*2/1,2 

*2/7,3 

• 

• 

• 

•*2/i.2/i 

ft/. 

*2/i  ,2/J+l 

*2/i  ,2/1+2 

b2n 

(2.27) 

Gauss  Elimination  is  then  used  to  solve  for  the  qt  in  terms  of  y(l),  y(2)  and  the 
constant  term.  This  gives 

qj  =  b\jy{\)  +  b2j  y(2)  +  bl}  j  =  1,2 ...  2n  (2.28) 

Equation  2.28  is  then  substituted  into  the  Kutta  condition  at  the  two  trailing  edges  to 
form  two  linear  simultaneous  equation  with  two  unknowns  y(l)andy(2)  .  Gauss 
Eimination  is  again  used  to  solve  for  the  Vorticity  distribution  y(l)andy(2)  and  these 
are  then  back  substituted  to  solve  for  the  qr 

3.  Computation  of  Velocity  and  Pressure  Distribution 

Once  the  q \  (J  =  1,2 ...  2n)  and  y ,(/*=  1,2)  are  solved,  the  velocities  at  all  the  panel 
control  points  can  be  easily  obtained.  The  normal  velocity  is  given  by 
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n  2n  n  In 

+  y(2)  ]T  Bny  +  V"  Sin| «(/)  -  0,  ]  (2.29) 

/■I  jmn+\  Jm  1  jmn+ 1 

while  the  tangential  velocity  is  given  by 

n  2n  n  2n 

(n  -  Yja,m+  Yj  + y(  + y(2)  X  ^ +  (2-3°) 

7-1  Jmn+l  Jm  l  7— n+1 

with 

1-1,2.../!  ;/-l, 

=  n  +  1  ...  2«  ;  =2 

The  normal  velocity  will  obviously  satisfy  equation  2.9  (used  to  check  the  code)  showing 
that  the  flow  tangency  condition  is  satisfied  while  the  the  total  velocity  will  be  given  by 
the  tangential  velocity. 

Vt0M  -  (V),,  /  =  1,2  ...  2n  (2.31) 

w’here 


<n  - 


w  n  2« 

(V1),  «  -  -  ]T  B^qj  +  y(l)^4"  +  >’(2)  X  4;  +  F^Cosl «(/)  -  6t ]  (2.32) 

7«i  jm, i+i  7-1  7-fl+i 

Substituting  equation  2.31  into  2.8  for  C,  with  (F'),  defined  as  in  equation  2.32,  the 
pressure  coefficient  at  the  /*  control  point  is 

(C,),  *  1  —  ( V'jj  ,  /—  1,2.. .2/i  (2.33) 

4.  Computation  of  Forces  and  Moments 

The  two-dimensional  aerodynamic  lifl  (Q,  drag  (CJ  and  pitching  moment  (C„) 
are  calculated  with  respect  to  the  global  frame  of  reference.  The  moment  coefficients 
are  computed  with  respect  to  the  respective  leading  edges.  For  the  code,  we  have 

n 

m  -  £(Cp)M+,  -  A',)  (2.34) 

/-I 
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'tit-ihimg&t 


CM-  -]T(C, 

Ml 

n 

cm  -  £(<:,),{  (A’h.,  -  XDXm, +{YW-  rjy«i, ) 

tel 

where  /  -  1,2  and  A„  A*,,  y„  Am,,  Km,  are  defined  as  before. 
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TiiatAt  lua 


(2.35) 

(2.36) 


III.  UNSTEADY  FLOW  PROBLEM  FORMULATION 


A.  GENERAL 

The  modification  work  for  the  unsteady  flow  is  much  more  involved.  It  requires  the 
following: 

1.  The  establishment  of  five  frames  of  reference  viz  one  fixed  inertia  frame  of  reference 
(global),  two  moving  local  frames  of  reference*  and  two  frozen  local  frames  of  ref¬ 
erenced 

2.  Reformulation  of  the  two  Kutta  conditions  which  are  coupled  non-linearly.  The 
solution  requires  an  iterative  procedure  to  compute  for  the  two  y(f).  There  are  two 
possible  solutions  to  the  Kutta  condition  due  to  the  quadratic  nature  of  the 
equations.  The  solution  which  ensures  that  the  product  of  the  tangential  velocities 
of  the  first  and  last  panels  of  each  airfoil  is  negative  is  accepted  as  the  solution. 

3.  The  creation  of  a  new  subroutine  (SUBROUTINE  NEWPOS)  which  transforms 
all  coordinates  in  either  of  the  two  respective  local  frames  of  reference  to  the  global 
frame  of  reference.  This  simplifies  the  definition  of  the  airfoil,  wake  element  and 
core  vortices  relative  global  geometries.  The  code  requires  that  the  airfoil  be  de¬ 
fined  with  respect  to  the  respective  frozen  local  frames  of  reference  once  only. 
Subsequent  time  dependent  airfoil  motion,  wake  panel  behaviour  and  core  vortices 
convection  are  computer  generated. 

4.  The  introduction  of  a  more  accurate  method  to  obtain  the  velocity  potential  by- 
integrating  the  velocity  over  smaller  panels  on  the  airfoil  without  having  to  store 
large  arrays  of  influence  coefficients  which  are  not  needed  for  satisfying  the  flow 
tangency  condition. 

5.  Extension  of  the  influence  coefficient  to  include  the  effects  of  the  second  airfoil  with 
its  own  peculiar  wake.  This  also  requires  an  introduction  of  an  additional  influence 
coefficient,  that  on  the  wake  element  due  to  the  wake  element  from  the  other 
airfoil. 

B.  FRAMES  OF  REFERENCE 

The  inertia  (global)  and  the  two  frozen  local  frames  of  reference  at  a  specified  time 
r*  are  defined  as  for  the  steady  flow  case.  The  twro  moving  local  frames  of  reference  have 
their  x-y  axes  as  for  the  frozen  local  frame  of  reference  but  this  frame  is  moving  with  the 
airfoil. 


4  Moving  local  frames  of  reference  are  used  to  satisfy  the  flow  tangency  equation  which  am¬ 
plifies  to  equation  2.9. 

5  Frozen  local  frames  of  reference  are  inertia  frames  (used  in  the  steady  flow  case)  of  reference 
used  to  convert  the  core  vortices  of  the  previous  time  steps  with  respect  to  that  time  step  frame  of 
reference  and  subsequently  transformed  to  the  current  time  step  frozen  frame  of  reference. 
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C.  UNSTEADY  FLOW  MODEL 
1.  Rigid  Body  Motion 

The  rigid  body  motion  for  the  two  airfoil  system  is  an  extension  to  the  single 
airfoil  system.  Both  airfoils  are  considered  to  have  a  mean  velocity  of  -  Vm,  with  a  time 
dependent  translational  velocity  of-|i/(/y  +  V([)J)  and  a  rotational  velocity  -0(4 
which  can  be  in  phase  or  out  of  phase;  /  andj  are  unit  vectors  in  the  respective  local 
frames  of  reference  and  O  is  positive  in  the  clockwise  direction.  The  flow  will  be  deter¬ 
mined  with  respect  to  the  moving  local  frame  of  reference.  The  flow  tangency  condition 
takes  its  simplest  form  in  the  moving  frame  of  reference.  The  flow  tangency  condition 
seen  from  this  frame  of  reference  will  satisfy  equation  2.9.  The  unsteady  stream  velocity, 
Vmm  is  made  up  by  the  vector  sum  of  a  mean  velocity  V„,  a  time  dependent 
translational  velocity  [(/(/)*'  +  V{f)J  ]  and  a  rotational  velocity  Q(/),  (Figure  4). 

vsneam  -  VJ.  Cos[  *(f)k  jT  +  Sin[  «(/)*  )j}  +  \  U(l)k7  +  V(l)j ]  +  Q {l)k  (y7  -  xj  )  (3.1) 

The  disturbance  potential,  of  necessity,  is  defined  with  respect  to  the  inertia  frame  of 
reference,  for  it  is  only  in  this  frame  that  the  flow  is  irrotational.  It  is  then  transformed 
to  the  moving  frame  of  reference  for  ease  in  treating  the  flow  tangency  condition.  The 
disturbance  potential  is  redcfmed  to  include  the  contributions  from  the  wake  panel  and 
the  core  vortices  from  both  airfoils. 

4>**4>s  +  <t>v  +  <i>w+<t>  cv  (3.2) 

The  total  velocity  is  then  written  : 


^  total  m  l  stream  "F  W*  (3.3) 

The  unsteady  Bernoulli's  equation  for  the  pressure  coefficients  on  the  airfoil  surface  is 
written  with  respect  to  the  moving  frame  of  reference.  Giesing  showed  this  to  be  written, 
in  our  notation  as : 


/  ^ stream  .2  .  ^ total  .2  _  2  3<t> 

'  Voo  v~  }  -  Vl  3' 

oo 


(3.4) 


where  Vtmm  and  Vn„,  are  defined  according  to  equations  3.1  and  3.3  respectively. 
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Figure  4.  Frames  of  Reference  for  Airfoil  1 
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2.  Wake  Model 

Recall  that  the  unsteady  flow  model  requires  an  additional  model  to  simulate 
the  continuous  shedding  of  vorticity  into  the  trailing  wake.  The  treatment  of  this  vortex 
shedding  process  follows  the  approach  of  Basu  and  Hancock.  The  shed  vorticity  takes 
place  through  a  small  straight  line  wake  element  attached  as  an  additional  panel  to  the 
trailing  edge  of  each  airfoil  with  a  uniform  vorticity  distribution  (y„(/)]*.  This  shed 
vorticity  panel  will  be  established  if  its  length  A (/)*  and  inclination  0 (/)*  to  the  respective 
local  frames  of  reference  satisfy  the  Helmholtz  theorem : 

A(0*ly*X0J*  +  r^0 -!■*_,(/)  (3.5) 

or  A  (/)  *  [  yjt)  ]*  -  r*_,(/)  -  rk(f)  -  SS(0  [  ?(/)*_,  -  y(/)* )  (3.6) 

where  SS  is  the  perimeter  of  the  airfoil  and  T*_,  and  y*_,  are  respectively  the  total  circu¬ 
lation  and  vorticity  strength  which  are  determined  at  the  previous  time  step 

At  the  next  time  step  r*l„  the  shed  vorticity  panel  will  be  detached  from  the 
trailing  edge  and  will  convect  downstream  as  a  concentrated  free  vortex  with  circulation 
A(0*  ly.(01*  or  r*_,(0  -  rt(0  at  the  resultant  local  velocity  of  the  fluid  particle.  This  wake 
convection  process  is  illustrated  in  Figure  5  where  the  airfoils'  subscripts  are  dropped 
without  loss  of  generality. 

In  the  code,  the  convection  of  the  core  vortices  is  broken  into  three  steps: 

1.  The  core  vortices  are  first  convected  using  the  resultant  absolute  velocity  with  re¬ 
spect  to  the  frozen  local  frame  of  reference  of  the  previous  time  step. 

2.  This  is  followed  by  a  transformation  to  the  current  frozen  local  frame  of  reference. 

3.  Finallv  this  is  transformed  to  the  global  frame  of  reference  by  using  the  new  sub¬ 
routine  (SUBROUTINE  NEWPOS). 

3.  Additional  Boundary’  Conditions 

The  unsteady  flow  model  has  now  introduced  an  additional  boundary  condition 
viz.  the  conservation  of  vorticity  (equation  3.5  or  3.6)  through  the  modeling  of  the  wake. 
However,  the  introduction  of  the  wake  creates  three  additional  unknowns  for  each 
airfoil,  that  is,  y „(/)».  A(/)»  and  ©(/)*.  As  such  two  additional  conditions  are  required  for 
each  airfoil  in  order  to  solve  the  system.  The  approach  suggested  by  Basu  and  Hancock 
is  extended  to  the  two  airfoil  case. 

1.  The  wake  panel  is  oriented  in  the  direction  of  the  local  resultant  velocity  at  the 
panel  midpoint. 


m 


life  l 


P 
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Vortex  Shedding  at  Time  Step  t^ 


Helmholtz's  theorem 
^k  +  **k  “  *"k-l 


^k*3  _ 

^k-4*“*"k-3 


(Vw)k 


tan0k  -  l-r\ 

(Lw)k 


^-<tk-Wv'l(l'w)k2  +  <Vw>k2J 


Figure  5.  Panel  Methods  Representation  for  Unsteady  Flow 
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tan  ©(/)* 


(3.7) 


\vlmk 

\vxMk 

where  [ ^ »(/))*  and  [Vi{t)\k  are  the  x  and  y  velocity  components  at  the  midpoint  of 
the  wake  panel  with  respect  to  the  frozen  local  frame  of  reference*. 

2.  The  length  of  the  wake  panel  is  propotional  to  the  magnitude  of  the  local  resultant 
velocity  at  the  panel  midpoint  and  the  step  size  of  the  time  step. 

m,  -  <f»  -  krow’+mow’  (3-*) 

We  now  have  all  the  necessary  equations  to  solve  for  the  system. 

D.  INFLUENCE  COEFFICIENTS 

As  with  the  boundary  condition,  additional  influence  coeflicients  need  to  be  defined 
as  a  result  of  the  wake  model.  The  definitions  in  the  single  airfoil  case  are  thus  extended 
as  follows: 

1.  More  A's  and  B's  Influence  Coeflicients 

Defining  NP3  and  NP4  as  the  panel  number  for  the  wake  for  the  first  and  sec¬ 
ond  airfoil  respectively  and  h  as  a  arbitrary  core  vortex,  the  following  additional  influ¬ 
ence  coefficient  are  required. 

•  (#!.vra)*  :  normal  velocity  component  induced  at  the  t*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  tk. 

•  (#lvw )*  :  tangential  velocity  component  induced  at  the  f*  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  /*. 

•  (#>«)*  ’•  normal  velocity  component  induced  at  the  f*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  tk. 

•  (5;_vw)»  tangential  velocity  component  induced  at  the  i*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  t„. 

•  :  x  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib¬ 
ution  on  the,/1*  panel  at  time  r*. 

•  ('ff.-njt  :  y  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib¬ 
ution  on  the/*  panel  at  time  rk. 

•  Ml vHjh  ■  x  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib¬ 
ution  on  the/*  panel  at  time  tk. 


6  It  is  important  to  note  that  the  velocities  at  the  midpoint  of  the  wake  panels  do  not  include 
the  effect  due  to  the  vorticity  distributed  along  itself  but  it  does  include  the  contribution  due  to  the 
wake  panel  of  the  other  airfoil. 
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•  y  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  source  distrib¬ 
ution  on  the/*  panel  at  time  tk. 

:  x  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis¬ 
tribution  on  the/*  panel  at  time  tk. 

(fyfnj*  •  y  velocity  component  induced  at  the  first  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis¬ 
tribution  on  the/*  panel  at  time  tk. 

'  x  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis¬ 
tribution  on  the  J*  panel  at  time  /«. 

:  y  velocity  component  induced  at  the  second  airfoil  wake  panel  midpoint 
with  respect  to  the  frozen  local  frame  of  reference  by  unit  strength  vorticity  dis¬ 
tribution  on  the/*  panel  at  time  tk. 

(AkJ)k  :  x  velocity  component  induced  at  the  h **  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  source  distribution  on  the 
/*  panel  at  time  tk. 

(AfJk  :  y  velocity  component  induced  at  the  h*  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  source  distribution  on  the 
/*  panel  at  time  tk. 

(5;/*  :  x  velocity  component  induced  at  the  h *  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  vorticity  distribution  on  the 
/*  panel  at  time  t„. 

(5U>  :  y  velocity  component  induced  at  the  h*  core  vortex  with  respect  to  the 
frozen  local  frame  of  reference  by  unit  strength  vorticity  distribution  on  the 
/*  panel  at  time  tk. 


’  x  velocity  component 
vorticity  distribution  on  the  wake 

($Lvn)*  :  y  velocity  component 
vorticity  distribution  on  the  wake 

:  x  velocity  component 
vorticity  distribution  on  the  wake 

’•  y  velocity  component 
vorticity  distribution  on  the  wake 


induced  at  the  h*  core  vortex  by  unit  strength 
panel  of  the  first  airfoil  at  time  tk. 

induced  at  the  h*  core  vortex  by  unit  strength 
panel  of  the  first  airfoil  at  time  tk. 

induced  at  the  A*  core  vortex  by  unit  strength 
panel  of  the  second  airfoil  at  time  t„. 

induced  at  the  h*  core  vortex  by  unit  strength 
panel  of  the  second  airfoil  at  time  t„. 


2.  New  C's  Influence  Coefficient 

The  single  airfoil  definition  is  extended  to  include  the  effects  of  the  second 

airfoil. 

•  (Q»(0)*  •’  normal  velocity  component  induced  at  the  i*  panel  control  point  by 

unit  strength  m'*  core  vortex  at  time  tk. 


ms- 
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•  (CUO).  :  tangential  velocity  component  induced  at  the  /'**  panel  control  point 
by  unit  strength  m*  core  vortex  in  the  wake  of  the  A*  at  time  /*. 

•  (C’xnAQ)*  '•  x  velocity  component  with  respect  to  the  frozen  local  frame  of  refer¬ 
ence  induced  at  the  first  airfoil  wake  element  by  unit  strength  m *  core  vortex  in  the 
wake  of  the  l*  at  time  it. 

•  :  y  velocity  component  with  respect  to  the  frozen  local  frame  of  refer¬ 
ence  induced  at  the  first  airfoil  wake  element  by  unit  strength  m4*  core  vortex  in  the 
wake  of  the  l*  at  time  t„. 

•  (Cj.^/))*  :  x  velocity  component  with  respect  to  the  frozen  local  frame  of  refer¬ 
ence  induced  at  the  second  airfoil  wake  element  by  unit  strength  m*  core  vortex  in 
the  wake  of  the  I*  at  time  i„. 

•  (Q ■  y  velocity  component  with  respect  to  the  frozen  local  frame  of  refer¬ 
ence  induced  at  the  second  airfoil  wake  element  by  unit  strength  m*  core  vortex  in 
the  wake  of  the  t*  at  time  rt. 

•  (Qj»(/))*  :  x  velocity  component  with  respect  to  the  frozen  local  frame  of  reference 
induced  at  the  h"1  core  vortex  by  unit  strength  m*  core  vortex  in  the  w>ake  of  the  l* 
at  time  /*. 

•  (Q.m(0)*  :  y  velocity  component  with  respect  to  the  frozen  local  frame  of  reference 
induced  at  the  h*  core  vortex  by  unit  strength  m*  core  vortex  in  the  w*ake  of  the  t* 
at  time  it. 

3.  Computation  of  the  Time-Dependent  Influence  Coefficient 

As  for  the  single  airfoil  ,  the  influence  coefficients  of  the  wake  element 
( .( and  are  computed  as  for  the  airfoil  influence  coefficients 

fij  and  Bij  with  the  subscripts  NP3.NP4  replacing  j. 

The  x  and  y  velocity  components  for  the  respective  frozen  local  frames  of  ref¬ 
erence  are  obtained  indirectly  by  first  calculating  for  the  global  X  and  Y  velocity.  These 
global  velocities  are  then  transformed  to  the  frozen  frames  through  a  simple  relationship 
by  the  use  of  the  respective  angle  of  attack.  It  can  be  easily  shown  that  the  velocity  with 
respect  to  the  frozen  frame  of  reference  is  related  to  the  global  velocity  by  the  following 
relationship. 

(O*  -  cos  <*(/)]*  -[GP^  sin  «(/)]* 


(Ok  «  [GOin  «(0]*  +  IGYU  cos  «(/)]* 


(3.9) 


w'here  V  denotes  a  generic  induced  velocity  with  respect  to  the  frozen  frame  and  the 
precedent  G  denotes  global  velocities. 
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The  global  ( GAt/nJ)t,{ GA^f^)k, ...  (G2ijrM>,)*,(G.,4')4,(G/4|,)*,(GZ?,’)*  and 
(GBfj)„  are  computed  using  equations  2.18  through  2.21  with  6 ,  set  to  zero  and  subscript 
i  appropriately  replaced. 

The  global  C's  coefficients  will  be  computed  identically  as  for  the  single  airfoil 
case.  The  results  are  repeated  here  for  the  sake  of  completeness. 


(GC,m 


CosjWDk  ~  (6m)k) 

M'bJk 


(3.10) 


(GCUWk 


Sinmk-(6m)k) 

M'lmh 


where  : 

(O*  *  sj[(xm,  -  x„y  +  (ym,  - 
xm,  -  1/2  (jf,  +  *,_,) 
ym,  -  1/2  (y,  +>^i) 

xm  =  x  coordinate  of  m *  core  vortex  at  time  t„ 
ym  -*  y  coordinate  of  m*  core  vortex  at  time  ik 


et 


ian~'( 


y~>-y< , 

*,*.  -  X,  > 


(<u 


tarr'{ 


ym,-ym 

xm,  —  xm  1* 


(3.11) 


Also  (GCjrpj^/))*,  (GCfo^,^/))*,  (GCfo*^/))*,  (GQrw>m(/))*,  (GCj^X/))*  and  (GC{^,(f))t  are 
computed  with  equation  3.10  with  6,  set  to  zero  and  the  subscript  i  appropriately  re¬ 
placed. 


E.  NUMERICAL  SOLUTION  SCHEME 
1.  The  Flow  Tangency  Condition 

The  unsteady  flow  tangency  equation  for  the  two  airfoil  system  is  a  simple  ex 
tension  of  the  single  airfoil  case. 


n  in  n  in 

+  !(?.»»»)/•  "J* 


1  J—n+) 


Jm  1  Jm, |+1 


+  lVw(l)]*(fi"v«)ft  +  tVw(2)l*(^W)ft  + 
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*-l  Jt-I 

|2jl(C(0)(r„.,(/)  -  rjO)l|j.i  +  &(c£,»)(r„-,(/)  -  r„(0)]Jw 


where  i*  1,2  ...  2n 


(3.12) 

where  i- 1,2  ...  2n  and  Vtlrmm  is  evaluated  by  equation  3.1  at  the  /•*  panel  control  point. 

Rearranging  the  equation  we  can  express  the  source  strengths  explicitly  as  a 
function  of  the  two  vorticity  strengths  and  a  constant  term  as  in  equation  2.28  where  : 


n  7n 

L.H.S.  - 


Jmn+\ 


First R.H.S.  - 


-  X/s;>*} 

j-\  ) 


(3.13) 


(3.14) 


Second  R.H.S.  =* 

Third  R.H.S.  « 


2  n 


(  Jmn+l  J 

- i’t)k  - 


(3.15) 


A-l  *-l 

-  jX)(C(0)(r„-,«  -  r„w)ij,.,  -  |gl(C»xr^,»  -  r„»)ijw  (3.16) 


where  i  *  1,2  ...  2n 

For  the  code,  the  above  is  implemented  in  a  2Ar  x  (2Ar  +  3)  matrix  in  SUB¬ 
ROUTINE  COEF  and  is  subsequently  solved  by  GAUSS  elimination  to  obtain  an  ex¬ 
plicit  expression  for  the  source  strength  in  terms  of  the  vorticity  distribution  and  a 
constant  term  as  in  the  steady  flow  case  as  follows : 
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qj  «  blj  •  y(l)  +  b2j .  y( 2)  +  b3j  ,  j  =  1,2 ...  2« 
where 


(3.17) 


61;  :  part  of  the  source  strength  which  depends  on  circulation  y(l)*. 
blj  :  part  of  the  source  strength  which  depends  on  circulation  y(2)*. 
b3j  :  part  of  the  source  strength  which  is  independent  of  the  circulation. 


Equation  3.17  is  critical  in  the  simplification  of  the  treatment  for  the 
non-linear  coupled  Kutta  condition. 

2.  The  Kutta  Condition 

For  the  unsteady  flow,  the  Kutta  condition  can  be  written  as  the  following  : 
For  the  first  airfoil 


»,)U  -  2( 


Em 

ot 


)k 


*=  255(1) 


y*(i)-y*-i(D 

Or—  i 


(3.18) 


For  the  second  airfoil 

(O l  -  (K*)l  -  2(  £  (*J»  -  i»  -  2(  )„ 


Si 


=  255(2) 


^(2)-»-i(2) 

lk  ~  lk- 1 


The  tangential  velocity  for  the  first  panel  can  be  written  thus: 


(3.19) 


2  n  n  2  n 

y;  _  -YJ(B'A)*  +  *'EA"l  +  +  A2)tYiA'i  + 


y- 1 


Jmtt+l 


^r.»«[y*-i(D  -  y*0)j  +  -  y,(2)l  +  •  ?1»  + 

A-l  k-l 

|El(c''-«Kr„-,(0  -  rm(0)lj,„  + 1  V«c;m(0xr^,(0  -  r„</))l}«  p.») 


Substituting  equation  3.17  into  3.20,  the  following  simplified  relation  for  V[  is  obtained. 
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V[  -  Dlxy{\)k  +  D2xy(2)k  +  D3X  (3.21) 

where 
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(3.22) 

Similar  expressions  can  be  obtained  for  the  other  trailing  edge  panels  as  Follows  : 

K  -  />1«V(1)*  +  D2ny{2)k  +  D\  (3.23) 

C,  -  ^ln+1v(l)*  +  Z>2„+ly(2)*  +  D\+x  (3.24) 

^  -  ^’(D*  +  i>22ny(2)A  +  Z?32n  (3.25) 

where  the  coefficients  Dl,  D2,  D3  are  obtained  as  per  equation  3.22.  We  require  the 
square  of  the  tangential  velocity  for  the  Kutta  condition 

(V[)1~{D\xy{\)k  +  D2xy(2)k  +  D3x)1 

~Dlly(l)l  +  D2]Y(2)l  +  D3l  (3.26) 

+  2D\xD2xy{i)ky{2)k  +  201,03, y(l)*  +  202,03,y(2)* 

(K')2«(Dlny(l)*+D2„y(2)*  +  Z)3„)2 

«012y(l)2  +  022y(2)2  +  032  (3.27) 

+  2D\„D2„y(l)ky(2)k  +  2Dl„D3„y(l)*  +  2D2„D3ny{2)k 
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For  the  first  airfoil  equations  3.26  and  3.27  are  then  substituted  into 
give,  for  the  code: 

AAA(\)y{l)2k  +  BBB(.\)y{2)\  +  CCC(l)y(l)*  + 

Z)Z)Z)(l)y(2)k  +  £££(l)y(l)Jty(2)*  +  FFF[l)  -  0 
where 

AAA(  1)  =  Dl]  —  Dl* 

BBB(  1)  *  D22-D22 

CCC(  1)  =  2[/)l,.Z?3»-Z?lw.Z?3w-  tS5!-  l 

lk  lk- 1 

DDD(l)  -  2{Z?2, .  Z?3,  -  D2n .  Z>3„) 

£££(1)  -  2(01 1 •  D2\  —  D\n*D2n\ 

£££(1)  -  Z>3?  -  D32n  +  2(  ) 

•k  ~  lk- 1 

A  similar  set  of  equation  can  be  obtained  for  the  second  airfoil  as  follows: 
AAA(  2)y(l)J  +  BBB{  2)y(2)|  +  CCC(2)y(l)k  + 

DDD(2)y(2)k  +  EEE(2)y(l)ky(2)k  +  fff(2)  -  0 
where 

AAA(2)  -  D\2n+i-D\\n 
BBBi 2)  «  D22n+,  -  D2\n 
CCC( 2)  -  2[/>l„+I .  D3„+1  -  D\ln .  £>3^] 

DDD(2)  -  2[f)2n+1 .  Z>3„+1  -  D2ln .  f>32„  -  ] 

£££(2)  -  2[Z)1W+, .  D2n+l  -  D\2n .  Z>22„J 


3.18  to 


(3.28) 


(3.29) 


(3.30) 
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rm 2)  -  />3«,  -  D>1  +  2(-^^f2)-)  (3-31) 

The  above  Kutta  equations  are  non-linear  and  coupled.  To  linearize  it,  we  adopt 
Newton's  method,  in  which  the  initial  iteration  is  linearized  about  the  values  of  the 
previous  time  step: 

vOS  -  rtOnl  +  Mfi-MlS)1  -  IKDl-.l1  + 

y(2)J  -  7(2);:!  +  WDl  -  (7(2>;i!  -  W2);.,l:  +  2y(2);;|«y(2);  (3.32) 

where: 

WDtevar*-'  ;  '/(!)*“  y(D*-i 
«5y(2)^y(2)r‘  ;  >’(2)* «  v(2)k_1 

with  k  denoting  the  iteration  counter.  After  dropping  quadratic  terms  of  <5y  and  col¬ 
lecting  like  terms  the  Kutta  condition  for  the  first  airfoil  simplifies  to  the  following: 

(2AAA(  l)y(I)*_,  +  CCC(  1)  +  £££(l)y(2)^1}«5y(i)Jt  + 

{2BBB(l)y(2)k_l  +  DDD(l)  +  £E£(l)y(l)k,1}3y(2)k  + 

■  AAA(l)y(l)l,  +  BBB(l)y(2)l,  +  CCC(l)y(l)*_,  + 

Z>Z>Z>(l)y(2)Jk_1  +  £££(  1 )'/( 1  )*_ i y(2)*_]  +  FFfl  1)  -  o 

(3.33) 

Similarly  the  linearised  Kutta  condition  for  the  second  airfoil  can  be  written: 
{2AAA(2)y(l)k_,  +  CCC(2)  +  £££(2)y(2),_,}<5y(l)*  + 

{2BBB(2)y(2)k.i  +  DDD(2)  +  £££(2)y(l)*.1}<5y(2)*  + 

AAA(2)y(l)2k.l  +  BBB(2)y(2)i_i  +  CCC(  2)y(l)*_,  + 

DDD( 2)y(2)Jt_1  +  £££(2)y(l)*.1y(2)*.l  +  FFF{2)  -  0 

(3.34) 
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The  above  two  linear  equations  with  two  unknowns  6y(l)k  and  Sy(2)k  can  now  be  easily 
solved  with  the  same  Gauss  routine.  The  results  are  then  back  substituted  into 
equations  3.32.  This  procedure  is  repeated  until  the  corrections  <5y(l)j  and  <5y(2);  are  less 
than  a  prescribed  tolerance. 

3.  Convection  of  the  Wake 

The  resultant  velocities  of  all  core  vortices  are  calculated  using  the  frozen  local 
frame  of  reference.  This  is  resolved  into  components  at  the  h*  core  vortex  as  follows: 

n  2n  n  2n 

<t4).  -  + Y  iA^)]t + ^Y^ +v(2,*X  + ,( • 7  ■  i* 

>1  /■«+ 1  >1  Jmn+ 1 

+  +  ly  wf  2)]*(  A'/m)* 

+  jXl(C^0Wrm-,(0  -  r«(0)]J/-t  +  jSKCSUOWr—iW  -  rM(/»]  j,.2 

m+h  m+h 

(3.35) 

n  In  n  2n 

(’'.)» -  + y  w1* + + Mi‘Y(-B>v)k + ^  J]k 

Jml  >11+1  >1  Jmn+ 1 

+  (yw(*)lk(^A,,Vw)fc  +  (Vw(2)I*(#ft,,VM)* 

+ 1  YjlicUOWm-At)  -  rm(/))l  L,  +  \  SKCjU0)a(r»-i(0  -  rmm  L2 

m+h  m+h 

(3.36) 


The  location  of  the  core  vortices  at  the  new  time  step  is  then  computed  with  respect  to 
the  current  time  step  frozen  local  frame  of  reference  and  then  transformed  to  the  new 
time  step  frozen  local  frame  of  reference.  These  coordinates  are  subsequently  trans¬ 
formed  to  global  coordinates  so  as  to  facilitate  the  calculation  of  the  influence  coeffi¬ 
cient. 

4.  Disturbance  Potential  and  Pressure  Distribution 

The  concept  of  disturbance  potential  introduced  in  the  steady  flow  has  played 
a  vital  role  in  that  it  facilitates  the  synthesis  of  the  flow  field  from  simple  flow  field  by 
simple  superposition  of  the  various  singularities  contributions.  However,  there  has  not 
been  a  requirement  to  solve  directly  for  the  disturbance  potential,  as  our  interest  was 
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on  the  disturbance  induced  velocity  which  is  the  spatial  derivative  of  the  disturbance 
potential.  In  our  solutions,  the  concept  of  influence  coefficient  has  allowed  a  direct 
evaluation  of  the  disturbance  velocity  thus  nullifying  the  requirements  of  obtaining  the 
disturbance  potential. 

The  treatment  for  the  unsteady  flow,  though  it  follows  the  same  procedure  as 
the  steady  flow'  case,  viz  the  influence  coefficient  to  obtain  the  disturbance  velocity,  still 
requires  the  disturbance  potential  for  the  computation  of  the  pressure  coefficient  as  can 
be  seen  in  equation  3.4.  Rewriting,  we  have: 


w* 


K- 


stream 


■Xfi-K' 


—  )& 


v. 


2  (#/)*  ~  (^<)fc-i 

V2  *k  -  lk-l 

oo 


/  «  1,2 ...  Zw 


(3.37) 


where  V„,„m  and  V„„,  are  defined  according  to  equations  3.1  and  3.3  respectively. 

From  equation  3.37  w*e  have  written  the  rate  of  change  of  <f>  by  a  backward  fi¬ 
nite  difference  approximation.  This  simplifies  our  iteration  procedures  tremendously  as 
the  <{>  from  the  previous  time  step  exists  at  the  current  computation.  The  computation 
of  the  disturbance  potenial  </>  is  obtained  through  tw'o  steps.  The  difference  in  potential 
<f)  from  upstream  at  infinity  to  the  leading  edge  is  computed  and  is  then  combined  with 
the  difference  in  the  potential  from  the  leading  edge  to  the  panel's  of  interest  control 
point. 

The  present  approach  differs  from  the  original  single  airfoil  approach  in  that  the 
velocity  potential  along  the  airfoil  surface  is  computed  via  a  finer  grid  using  Gaussian 
quadrature?.  This  modification  improves  the  results*  but  the  cost  of  doing  it  is  a  longer 
computation  time.  The  definition  of  infinity  has  also  been  set  to  100  chord  lengths  in 
comparison  to  10  chord  lengths  used  in  the  original  code.  However,  the  number  of 
computation  points  and  the  first  panel  length  from  the  leading  edge  upstream  have  not 
been  changed. 

For  completeness,  the  computation  of  the  disturbance  potential  from  the  lead¬ 
ing  edge  to  infinity  is  included  here  for  both  airfoils.  It  is  essentially  the  same  as  the 
single  airfoil  case  except  that  the  disturbance  potential  is  now  computed  relative  to  the 
global  frame  of  reference  instead  of  the  airfoil  fixed  frame  of  reference  as  used  in  the 

7  Each  panel  is  subdivided  into  4  additional  sub-panels  and  the  tangential  velocity  is  computed 
and  integated  over  these  smaller  panels  with  a  weighting  function  to  get  the  disturbance  potential. 

8  One  can  check  the  improvement  by  computing  the  y(r)  obtained  through  the  difference  be¬ 
tween  the  trailing  edges  and  compares  them  with  that  obtained  through  the  Kutta  condition. 


33 


original  code.  We  begin  by  selecting  a  straight  line  extending  upstream  in  the  direction 
parallel  to  V„.  The  length  of  the  line  is  set  at  100  chord  lengths.  This  line  is  divided  into 
z  panels  with  the  first  element  at  the  leading  edge  set  equal  to  the  single  airfoil  case  for 
the  purpose  of  ensuring  that  the  panel  size  is  comparable  to  the  airfoil  panel  size.  The 
panel  size  is  then  subsequently  increased  to  take  advantage  of  the  inversely  decaying 
induced  velocities  at  larger  distances.  Using  subscript  f  to  denote  these  panel  mid-points, 
we  define  the  following  influence  coefficients: 

•  (A^)„  :  normal  velocity  component  induced  at  the/*  panel  control  point  by  unit 
strength  source  distribution  on  the  /panel  at  time  r*. 

•  (A'fjt  :  tangential  velocity'  component  induced  at  the/*  panel  control  point  by  unit 
strength  source  distribution  on  the  /panel  at  time  r*. 

•  :  normal  velocity  component  induced  at  the  /  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  /panel  at  time  r*. 

•  (££,)*  :  tangential  velocity  component  at  the/  panel  control  point  by  unit  strength 
vorticity  distribution  on  the  /panel  at  time  tk. 

•  (fyvn)*  ■  normal  velocity  component  induced  at  the/*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  r*. 

•  (tyvn),  :  tangential  velocity  component  induced  at  the/  panel  control  point  by 
unit  strength  vorticity  distribution  on  the  wake  panel  of  the  first  airfoil  at  time  r*. 

•  ify v»)*  :  normal  velocity  component  induced  at  the/  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  r*. 

•  (®Lvw)»  tangential  velocity  component  induced  at  the/*  panel  control  point  by  unit 
strength  vorticity  distribution  on  the  wake  panel  of  the  second  airfoil  at  time  r*. 

•  (Q«.(0)>  •'  normal  velocity  component  induced  at  the/ panel  control  point  of  the 
/'*  by  unit  strength  m*  core  vortex  at  time  r*. 

•  (Q*(0)*  •'  tangential  velocity  component  induced  at  the  /  panel  control  point 

of  the  /'*  by  unit  strength  m,k  core  vortex  at  time  r*. 

The  tangential  velocity  at  the/  panel  written  as  for  the  code  is  then: 


>1  >1  >*+1 

fcUUD  -  r*0)l + -  »(2)!  + 

k~\  k- 1 


(3.38) 


valid  for  f»  1,2  ...  z.  The  disturbance  potential  at  the  airfoil  leading  edge  is  the  sum  of 
the  products  of  the  disturbance  induced  velocity  at  each  panel  and  the  panel  length  . 


<*  J.  -  -  SC/>«  •  (•‘■j-  -  A»1  (i-39j 

/-I 

The  integral  over  the  airfoil  surface,  as  stated  before,  is  now  done  over  a  finer 
grid.  This  requires  the  computation  of  the  disturbance  induced  velocity  over  the  smaller 
grids  within  each  panels.  Defining  the  total  finer  grid  points  in  one  panel  by  P9,  we  de¬ 
fine  first  the  refined  influence  coefficient 

MU 

tangential  velocity  induced  at  the  i*  panel  p*  node  due  to  unit  strength 
source  distribution  on  the  j*  panel  at  time  tk 

The  other  refined  influence  coefficients  have  the  same  definition.  The  tangential  com¬ 
ponent  of  the  disturbance  induced  velocity  at  the  i*  panel  m1*  node  is  then 

In  n  2  <i 

to*  -  -  yw?A + HD»y 4. + v(2).y,; + 

l  >»l  /■»+ 1 

(1)  -  ».(')!  +  -  y,(2)l  + 


*-l  *-l 

|Xj(«i,«)(r„.l(o  -  r„(/»i  jtal  +  |2[(ci,,(/)Krm.,(0  -  r„(0)l|M  (3.40) 

which  is  valid  for  i«  1,2  ...  2n  ;  p«  1,2  ...  P.  Performing  the  line  integration  by  summa¬ 
tion,  the  disturbance  potential  at  the  /*  nodal  point  is  then 


(<t>nodtl)k-(K)k+/  jj(V*)krJJ+iWp),  for  n>i^ik 


(3.41) 


~  7  '  forn>i*ik 

3^-1 


9  This  includes  the  end  points. 
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where  rJthl  denotes  the  panel  length  and  IV,  the  Gaussian  quadrature  weighting  function. 

(3-42) 


’m,  ~  s !<*»,-*/ HYm-rf 

Finally  the  disturbance  potential  at  the  i *  panel  control  point  is : 

{4>dk  “  l/2(  {4>„cHiei)k  +  (4>nod*l+\)k]  .  >  *s  1.2  ...  2 n 


(3-43) 


5.  Computation  of  Forces  and  Moments 

The  C„  Cd  and  Cm  about  the  leading  edge  of  each  airfoil  are  calculated  in  exactly 
the  same  way  as  it  is  done  for  the  steady  flow  problem  by  integrating  the  pressure  dis¬ 
tribution  (See  section  D-4  of  Chapter  2). 
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IV.  DESCRIPTION  OF  COMPUTER  CODE  USPOTF2 


A.  PROGRAM  STRUCTURES  AND  CAPABILITIES 

1.  Restrictions  and  Limitations 

The  restrictions  and  limitations  listed  in  the  single  airfoil  documentation  [Ref. 
1]  still  apply  for  the  two  airfoil  case.  Some  efforts  were  made  to  optimise  storage  re¬ 
quirements  versus  computational  repetitions;  however  as  in  the  original  code  there  is  still 
much  room  for  improvement.  Krainer  [Ref.  6  ]  has  improved  the  original  code  by 
combining  subroutines  and  reducing  repeated  computation  by  introducing  additional 
common  variables.  Some  of  his  improvements  have  been  added  to  this  program. 

The  computer  system  used  is  the  Naval  Postgraduate  School  IBM  3033AP.  The 
current  program  Axes  the  maximum  number  of  airfoil  panels  to  20010  and  the  maximum 
allowable  time  steps  to  200.  The  computer  currently  has  to  be  run  with  a  minimum 
storage  requirement  of  2  Mbytes.  A  detailed  computing  time  study  for  the  program  is 
not  undertaken,  for  it  not  only  changes  with  the  number  of  nodal  points  selected  but 
also  with  the  time  step  increment  as  this  has  an  effect  on  the  rate  of  convergence  for  the 
wake  panel  iteration.  An  order  of  magnitude  is  given  for  one  case  run  to  give  an  ap¬ 
preciation  of  the  computing  time  required.  The  system  currently  requires  a  total  CPU 
run  time  of  200  seconds  using  an  optimising  compiler  for  a  step  input  with  a  0.025 
non-dimensionalised  time  increment  for  26  time  steps. 

At  present,  the  program  can  run  for  two  airfoils  set  at  arbitrary  distance  and  at 
different  angles  of  attack  undergoing  any  of  the  following  motions: 

1.  In-phase  and  out-of-phase  Step  Input 

2.  In-phase  and  out-of-phase  Modified  Ramp  Input 

3.  In-phase  and  out-of-phase  Translational  Harmonic  Oscillation 

4.  In-phase  and  out-of-phase  Rotational  Harmonic  Oscillation 

5.  Sharp  Edge  Gust  Field  Penetration 

2.  Current  Structures  of  USPOTF2  Main  Program 

The  flow  logic  for  USPOTF2  (Unsteady  Potential  Flow  for  2  airfoils)  is  illus¬ 
trated  in  Figure  6.  It  is  essentially  divided  into  three  major  modules,  namely,  the  input 

10  These  are  combined  total  panels  for  both  airfoils. 
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-  problem  setup,  Steady  flow  solution  with  its  associated  output  and  the  Unsteady  flow 
solution  with  its  associated  output. 

The  first  module  includes  subroutines  INDATA,  SETUP,  BODY  and  NACA-15. 
The  main  functions  of  this  module  are: 

1.  Set  up  the  problem  formulation  by  reading  in  the  necessary  values  and  flag  setting 
from  filecode  l. 

2.  If  the  flag  for  a  NACA  4-digit  or  5-digit  of  type  230XX  is  set,  this  module  calls 
subroutines  BODY  and  NACA45  to  compute  the  local  panel  coordinates. 

3.  If  the  flag  for  a  NACA  4-digit  or  5-digit  of  type  230XX  is  not  set,  this  module  then 
proceeds  to  read  in  the  local  panel  coordinates  from  filecode  1. 

4.  The  local  slopes  and  the  airfoil  perimeters  are  then  calculated  in  preparation  for  the 
next  module. 

The  second  module  calculates  the  Steady  flow  solution.  It  calls  subroutines 
NEWPOS,  INFL,  COEF,  GAUSS,  KUTTA,  VELDIS  and  FANDM.  The  main  func¬ 
tions  of  this  module  are: 

1.  Transform  the  local  airfoil  coordinates  into  global  coordinates  and  compute  the 
influence  coefficients. 

2.  Set  up  flow  tangency  equation  as  per  equation  2.27  and  solve  for  the  source 
strengths  (q)  as  a  function  of  the  vorticity  strengths  [y(/)J  and  a  constant  part  where 
j  “  1,2  ...  in  and  /**  1,2. 

3.  Set  up  and  solve  the  Kutta  condition  as  per  equations  2.25  and  2.26  for  the 
vorticity  strengths  and  back  substitute  to  get  the  source  strengths. 

4.  Compute  the  total  tangential  velocity  in  the  moving  frame  and  compute  the  dis¬ 
turbance  potential  (<£,)  and  pressure  coefficient  [(QJ  (i  *  1,2  ...  2n). 

5.  Finally,  compute  the  aerodynamic  coefficients  of  forces  and  moments  - 
C  i  Q  i  C„. 

6.  This  program  can  terminate  in  this  module  after  all  the  steady  flow  parameters  are 
obtained  without  necessarily  running  the  Unsteady  flow  code. 

The  third  module  calculates  the  Unsteady  flow  solution.  It  calls  subroutines 
NEWPOS,  INFL,  COEF,  GAUSS,  KUTTA,  TEWAK,  PRESS,  FANDM  and 
CORVOR.  The  main  functions  of  this  module  are: 

1.  For  the  particular  unsteady  motion,  to  compute  the  initial  time  step  and  new  airfoil 
orientation  with  its  associated  local  body  velocities. 

2.  Introduce  the  wake  panel  and  assume  an  initial  length  and  orientation  to  begin  it¬ 
erations. 

3.  Transform  all  local  coordinates  and  panel  slopes  to  global  coordinates  and  global 
panel  slopes. 
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4.  Update  influence  coefficient  and  set  up  flow  tangency  equation. 

5.  Solve  for  the  source  strengths  in  terms  of  the  vorticity  strengths  and  the  constant 
part. 

6.  Invoke  the  required  Kutta  condition  and  solve  for  the  vorticity  strengthsil;  back 
substitute  to  get  the  source  strengths. 

7.  Update  the  local  velocities  of  the  wake  element  and  compute  new  length  and  ori¬ 
entation  •••  check  for  convergence. 

8.  If  wake  element  velocities  have  not  converged,  iterate  with  the  new  wake  element 
geometries  until  convergence  12. 

9.  Compute  the  total  velocity,  disturbance  potential  and  the  pressure  coefficient. 

10.  Compute  the  aerodynamic  lift,  drag  and  moment  coefficient. 

11.  Adjust  time  step  either  through  an  interative  inputs  or  by  a  automatic  time  in¬ 
crement. 

12.  Compute  resultant  local  velocities  of  the  core  vortices. 

13.  Convect  the  core  vortices  by  the  procedure  described  in  Chapter  III  section  E-3 


11  There  are  two  possible  solutions;  only  the  vorticity  distribution  that  ensures  that  the  prod¬ 
uct  of  the  relative  tangential  velocities  of  the  upper  and  lower  panels  of  the  trailing  edge  is  negative 
is  accepted  as  solution. 

12  The  convergence  criterion  is  user  specified  through  input  data  for  TOL. 

13  This  is  accomplished  by  setting  TADJ  to  be  non-zero  and  is  intended  for  use  in  conjunction 
with  the  viscous  flow  program. 
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B.  DESCRIPTION  OF  SUBROUTINES 

1.  Subroutine  BODY 

This  subroutine  is  called  for  the  purpose  of  obtaining  the  local  (x,y)  coordinates 
of  either  a  NACA  XXXX  or  230XX  type  airfoil.  It  is  called  by  subroutine  SETUP  and 
it  in  turn  calls  subroutine  NACA45  to  obtain  the  airfoil  thickness  and  camber  distrib¬ 
utions. 

2.  Subroutine  COEF 

This  subroutine  was  originally  intended  for  the  unsteady  flow  calculation  for  the 
single  airfoil  case.  It  is  now  modified  to  include  the  steady  flow  calculations.  Its  pur¬ 
pose  is  to  set  up  the  flow  tangency  matrix  as  in  equation  2.27  for  the  steady  flow  and 
equations  3.13  through  3.16  for  the  unsteady  flow  case.  These  matrices  are  necessarily 
set  up  in  this  way  so  that  the  source  strengths  can  be  solved  in  terms  of  the  vorticity 
strengths  and  a  constant  by  subroutine  GAUSS  as  a  linear  system  with  three  right  hand 
sides.  It  is  called  by  the  MAIN  program. 

3.  Subroutine  COFISH  (deleted  for  the  two  airfoil  case) 

This  subroutine  was  originally  set  up  to  serve  the  same  function  as  subroutine 
COEF  for  the  steady  case.  It  is  now  deleted  for  computational  efficiency. 

4.  Subroutine  CORVOR 

This  subroutine  is  called  for  the  purpose  of  obtaining  the  convective  velocities 
for  all  the  wake  core  vortices  with  respect  to  the  frozen  local  frame  of  reference  of  the 
current  time  step  in  accordance  with  equations  3.35  and  3.36  where  all  the  influence 
coefficients  are  now  locally  computed  (previously  done  in  subroutine  INFU)  to  save 
storage  requirements  since  these  are  only  required  in  this  subroutine.  The  local  influence 
coefficients  are  obtained  indirectly  by  first  obtaining  the  global  influence  coefficients  and 
transforming  them  viz  equation  3.9.  This  subroutine  is  called  by  the  main  program 
nearing  the  end  of  the  unsteady  flow  calculations  before  starting  a  new  time  step. 

5.  Subroutine  FANDM 

This  subroutine  is  intended  to  calculate  the  overrail  lift  coefficient,  drag  coeffi¬ 
cient  and  moment  coefficient  about  the  leading  edge  for  the  two  airfoils.  In  order  to 
preserve  the  option  of  obtaining  the  x  and  y  forces  with  respect  to  the  respective  local 
frames  of  reference,  the  original  method  rather  than  equations  2.34  through  2.36  is  im¬ 
plemented.  This  subroutine  is  called  by  the  MAIN  program  in  both  the  steady  and  un¬ 
steady  flow  computations  immediately  after  computing  the  pressure  coefficient. 
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6.  Subroutine  GAUSS 

This  subroutine  is  the  standard  linear  system  solver  that  employs  the  well- 
known  Gaussian  elimination  with  partial  pivoting  and  operates  simultaneously  on  a  user 
specified  number  of  right-hand-sides.  It  is  called  by  the  MAIN  program  in  both  the 
steady  and  unsteady  flow'  calculations.  In  order  to  use  GAUSS,  the  coefficients  of  the 
augmented  matrix  must  be  set  up  so  that  GAUSS  will  return  the  solutions  replacing  the 
corresponding  columns  of  the  augmented  matrix  that  were  initially  occupied  by  the 
right-hand-sides.  The  coefficient  set-ups  are  done  by  subroutine  COEF  for  both  steady 
and  unsteady  flow  problems. 

7.  Subroutine  INDATA 

This  subroutine  is  intended  to  read  in  the  first  three  to  five  cards  of  the  input 
data  depending  on  whether  I  FLAG  «  0.  The  first  three  cards  contain  some  description 
of  the  airfoil  type,  problem  definition,  I  FLAG  information  as  well  as  the  number  of 
lower  and  upper  panels.  If  I  FLAG  ■  0,  it  will  treat  the  airfoils  as  NACA  type  airfoils 
and  will  proceed  to  read  the  NACA  number  and  to  calculate  the  thickness  parameters 
that  will  be  required  by  subroutine  NACA45.  This  is  the  first  subroutine  called  by  the 
MAIN  program. 

8.  Subroutine  INFL 

This  subroutine  generates  most  of  the  influence  coefficients  that  are  needed  and 
shared  by  the  different  subroutines.  It  has  been  modified  to  include  the  steady  flow  case 
for  the  purpose  of  reducing  repeated  computations.  It  utilises  the  known  relative  ge¬ 
ometrical  parameters  of  the  singularities  to  cany  out  computation  based  on  equations 
2.18  through  2.21  for  the  steady  flow'  calculations  and  including  3.9  through  3.11  for  the 
unsteady  flow  case.  The  MAIN  program  calls  this  subroutine  in  every  iteration  cycle 
of  each  time  step  so  that  the  time  dependent  influence  coefficients  can  be  updated  as 
and  when  neccessary.  Time  independent  coefficients  are  computed  once  in  the  entire 
flow'  solutionsH.There  are  also  time  dependent  influence  coefficients  that  are  independ¬ 
ent  of  the  iterative  cycle  to  obtain  the  wake  panel  orientation^.  Those  influence  coef¬ 
ficients  involving  the  wake  core  vortices  are  also  independent  of  the  iteration  cycle  and 
are  also  updated  once  in  each  time  step;  but  the  process  of  the  update  is  more  compli¬ 
cated.  This  is  due  to  the  process  that  the  wake  need  first  to  be  convected  with  respect 
to  the  previous  time  step  frozen  frame  of  reference,  transformed  to  the  present  frame  of 

14  These  are:  on  the  airfoil  panel  by  the  panels  on  the  same  airfoil. 

15  These  are:  on  the  airfoil  panel  by  the  panels  on  the  other  airfoil. 
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reference  and  finally  transformed  to  the  global  frame  of  reference  by  subroutine 
NEWPOS.  Finally  the  influence  coefficients  involving  the  wake  panels  are  calculated 
as  frequently  as  the  number  of  iterations  take  to  find  a  converged  solution. 

9.  Subroutine  KUTTA 

This  subroutine  is  intended  to  solve  the  Kutta  equation.  It  has  been  modified 
to  include  the  steady  flow  solution.  After  the  source  strengths  have  been  determined  by 
subroutines  COEF  and  GAUSS  in  terms  of  the  vorticity  strengths  and  a  constant,  this 
subroutine  invokes  the  Kutta  condition  as  in  equations  2.2S,  2.26  for  the  steady  flow 
case  and  3.33,  3.34  for  the  unsteady  flow  case.  The  two  linear  equations^  are  then 
solved  again  with  Gaussian  elimination  to  obtain  the  vorticity  distributions.  For  the 
unsteady  case,  we  add  in  the  additional  requirement  of  finding  the  product  of  the 
tangential  velocities  at  the  upper  and  lower  trailing  edges.  These  are  demanded  to  be 
negative  as  there  are  basically  two  possible  solutions  to  the  Kutta  equations  due  to  its 
original  quadratic  nature. 

10.  Subroutine  NACA45 

This  subroutine  is  intended  to  calculate  the  camber  and  thickness  distribution 
of  the  NACA  4-digit  and  the  XACA  5-digit  airfoils  of  type  230XX  which  share  common 
thickness  distributions  with  the  4-digit  airfoil  having  the  same  thickness  to  chord  ratio. 
This  subroutine  is  called  by  subroutine  BODY  which  is  in  turn  called  by  subroutine 
SETUP  and  in  turn  called  by  the  MAIN  program. 

11.  Subroutine  NEWPOS 

This  is  a  new  subroutine  introduced  as  a  result  of  setting  up  the  five  frames  of 
reference.  Most  of  the  coordinates  computation  is  done  with  respect  to  the  respective 
local  frame  of  reference.  Its  purpose  then  is  to  transform  all  coordinates  in  the  respec¬ 
tive  local  frames  of  reference  to  the  global  frame  of  reference^.  This  is  necessary  for  the 
computation  of  the  influence  coefficients  as  well  as  in  the  convection  of  the  wake  core 
vortices.  It  facilitates  the  simple  requirements  of  defining  the  airfoil  once  only  with  re¬ 
spect  to  its  local  frame  of  reference.  Orientation  and  displacement  of  the  airfoil  in  the 
two-dimensional  plane  is  henceforth  calculated  through  this  subroutine.  This  subroutine 

16  For  the  unsteady  case,  the  original  two  equations  are  non-linear  and  were  subsequently 
linearised  in  the  discussion  in  Chapter  3. 

n  It  has  a  secondary  function  to  obtain  the  slopes  for  the  airfoil  panels  and  the  wake  element 
panels. 
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is  called  from  the  MAIN  program  in  several  locations  immediately  after  the  local  coor¬ 
dinates  are  computed. 

12.  Subroutine  PRESS 

This  subroutine  calculates  the  pressure  distribution  over  the  airfoil  panels  after 

the  iterative  solution  for  the  unsteady  flow  problem  has  successfully  met  the  convergence 

criterion.  It  first  computes  the  tangential  velocities  at  all  panel  control  points  using 

P  +  1 

equation  3.40  (with  p  —  — - — ),  then  performs  the  disturbance  potential  evaluation  at 
the  current  time  step  according  to  equations  3.41  through  3.43.  Together  with  the  dis¬ 
turbance  potential  data  obtained  from  the  previous  time  step,  it  calculates  the  pressure 
distribution  using  equation  3.37. 

13.  Subroutine  SETUP 

This  subroutine  sets  up  the  local  panel  nodal  coordinates  for  MAIN  program 
by  reading  the  fourth  through  seventh  data  sets  of  the  input  file  if  IFLAG  *  1  is  set. 
It  skips  the  data  reading  if  IFLAG  =  0  and  proceeds  to  set  up  the  node  distribution  and 
calls  subroutine  BODY  to  calculate  the  airfoil  local  coordinates.  The  node  distribution 
adopts  a  cosine  formula  in  order  to  have  closely  packed  panels  toward  the  leading  and 
trailing  edges  for  improvements  in  solution  accuracy.  Regardless  of  how  the  nodal  co¬ 
ordinates  are  obtained  ,  subroutine  SETUP  determines  the  local  panel  slopes  and  airfoil 
perimeter  length. 

14.  Subroutine  TEWAK 

This  subroutine  is  intended  to  calculate  the  resultant  velocity  components  at  the 
mid-point  of  the  shed  vorticity  panel  with  respect  to  the  current  frozen  frame  of  refer¬ 
ence.  These  velocity  components  are  necessary  to  ensure  that  the  correct  shed  vorticity 
panels'  length  and  orientation  are  established.  This  is  the  governing  criterion  for  the  it¬ 
erative  solution  scheme  of  the  unsteady  flow  case.  This  subroutine  is  called  by  the 
MAIN  program  at  every  iteration  cycle  of  each  time  step  for  the  unsteady  flow  calcu¬ 
lation. 

15.  Subroutine  VELD1S 

This  subroutine  essentially  performs  the  same  function  as  subroutine  PRESS  for 
the  steady  flow  case.  It  is  thus  redundant  and  could  be  deleted  with  some  modification 
work  necessary  for  subroutine  PRESS.  This  work  has  been  implemented  by  Krainer  for 
the  single  airfoil  case.  This  subroutine  is  called  by  the  MAIN  program  in  the  steady  flow’ 
computation. 
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C.  INPUT  DATA  FOR  PROGRAM  USPOTF2 


The  Input  data  is  similar  to  that  for  the  single  airfoil  case.  Program  USPOTF2  reads 
its  input  data  from  filecode  1.  An  example  of  an  input  data  file  is  attached  in  Appendix 
B  for  the  case  when  the  airfoil  nodal  coordinates  are  input  by  the  user.  Computer  gen¬ 
erated  airfoil  coordinates  are  another  option  that  can  be  selected  if  the  airfoil  chosen 
belongs  to  the  family  of  NACA  4-digit  or  5-digit  airfoils  of  type  230XX.  To  do  this,  the 
I  FLAG  parameter  is  set  to  zero  in  the  first  item  of  the  4*  set  of  data  card  and  replace 
the  5“  and  6*  set  of  cards  by  two  cards  containing  the  particular  airfoil  NACA  number 
for  the  two  airfoils  using  format  (15).  Figure  7  contains  an  itemised  description  of  the 
sequential  input  variables. 

D.  OUTPUT  DATA  FROM  PROGRAM  USPOTF2 

The  Output  data  is  similar  to  that  for  the  single  airfoil  case.  Appendix  C  contains 
a  sample  output  data  generated  by  using  the  input  data  set  from  Appendix  B.  Due  to 
the  repetitive  nature  of  the  output  as  a  function  of  time,  only  selective  time  set  data  are 
shown.  The  output  data  file  begins  with  writing  out  what  the  program  has  read  from 
the  data  file  followed  by  the  computed  nodal  coordinates  only  if  they  are  generated  by 
the  program;  otherwise  it  proceeds  to  write  the  airfoil  computed  perimeter  length.  The 
next  set  of  output  data  are  the  steady  flow’  solution  parameters  of  distributed  source 
strengths,  vorticity  strengths,  pressure-velocity  distribution,  force-moment  coefficient, 
potential  at  the  control  nodal  points  and  the  potential  at  the  leading  edge.  In  addition 
some  output  is  given  to  allow  an  assessment  of  the  accuracy  of  the  flow'  solution.  This 
consists  of  the  normal  component  of  the  velocities  at  the  panel  control  points  and  the 
integral  of  the  disturbance  tangential  velocity  along  the  airfoil  contour.  If  only  the 
steady  flow'  solution  is  required,  the  output  terminates  here. 

For  unsteady  flow,  in  addition  to  the  above,  the  unsteady  flow  parameters  are  also 
printed.  This  includes,  at  every  time  step,  the  iterative  printout  for  the  convergence  of 
the  length  and  orientation  of  the  wake  element  panel,  the  unsteady  solution  parameters 
on  the  airfoil  similar  to  that  of  the  steady  flow  parameters  as  well  as  the  trailing  wake 
core  vortices  data.  Again  the  same  checking  mechanisms  are  inserted  on  the  airfoil 
panels.  A  detailed  explanation  on  the  output  variable  names  is  listed  in  Figure  8.  All 
output  data  are  non-dimensional  quantities. 


Data  Sat  #1 

Format  (IS)  -  1  data  card. 

ITITLE 

•  Number  of  title  cards  to  be  used  in  Data  Set  #2. 

Data  Sat  #2 

Format  ( 20A4)  -  ITITLE  data  cards. 

TITLE 

-  Headings  printed  on  output  for  case  run 

identification. 

Data  Set  #3 

Format  (I5.2F10. 6)  -  1  data  card. 

NAIRFO 

-  Number  of  airfoil,  in  this  case  *  2. 

XSHIFT 

-  Relative  X  distance  of  the  2  airfoil’s  pivot  position 
with  respect  to  the  airfoil  global  coordinate  system. 

YSHIFT 

-  Relative  Y  distance  of  the  2  airfoil's  pivot  position 
with  respect  to  the  airfoil  global  coordinate  system. 

Data  Sat  #4 

Format  (315)  -  1  data  card. 

IFLAG 

-  0  if  airfoil  is  NACA  XXXX  or  230XX. 

-  1  otherwise. 

NLOWER 

-  Number  of  panels  used  on  both  airfoil  lower  surface. 

NUPPER 

-  Number  of  panels  used  on  both  airfoil  upper  surface. 

Data  Set  #5 

Format  (6F10. 6)  if  IFLAG  »  1  -  variable  data  cards. 

x(I),y(I) 

-  local  non -dimens ionalised  x-nodal  followed  by  y-nodal 

coordinates  for  airfoil  1.  Total  of  Nlower+Nupper+1 

nodal  points  divided  into  6  points  per  data  card. 

Data  Sat  #6 

Format  (6F10. 6)  -  variable  data  cards 

*(i),y(D 

-  local  non-dimensionalised  x-nodal  followed  by  y-nodal 

coordinates  for  airfoil  2.  Total  of  Nlower+Nupper+1 

nodal  points  divided  into  6  points  per  data  card. 

Figure  7.  List  of  Input  Variables. 
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Data  Sat  #7  Format  (6F10. 6,13)  -  1  data  card 

ALPI(l)  -  Initial  angle  of  attack  (AOA)  for  airfoil  1  in  degree 

ALPI(2)  -  Initial  angle  of  attack  (AOA)  for  airfoil  2  in  degree. 

DALP  -  Absolute  change  in  AOA  in  degree  for  non  oscillatory 

motions. 

-  Maximum  amplitude  of  AOA  change  in  degree  for 
rotational  harmomic  motions. 

TCON  -  Non-dimensional  rise  time  (F^t/c)  of 

AOA  for  motion  involving  modified-ramp  change  in  AOA. 

FREQ  -  Non-dimensional  oscillation  (a>c/FQO) 

for  harmonic  motions. 

PIVOT  -  The  length  from  the  leading  edge  to  the  pivot  point 

for  the  local  system. 

IPKASE  -  Flag  for  in-phase  and  out-of-phase  motion. 

-  0  out-of-phase  motion. 

-  1  in-phase  motion. 

Data  Set  # 8  Format  (8F10. 6)  -  1  data  card. 

U6UST  -  Magnitude  of  non-dimensional  gust  velocity  along 

global  x-direction. 

VGUST  -  Magnitude  of  non-dimensional  gust  velocity  along 

global  y-direction. 

DELHX(l)  •  Non-dimensional  translational  chordwise  amplitude 

for  airfoil  1 

DELHX(2)  -  Non-dimensional  translational  chordwise  amplitude 
for  airfoil  2 


Figure  7  (cont'd) 
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Data  Set  #8 
DELHY(l) 

DELHYC2) 

PHASE(l) 

PHASEC2) 

Data  Set  if 9 
TF 

DTS 

TOL 

TADJ 

SCLA 

N6IES 


Figure  7  (coat'd) 


(Cont'd) 

*  Non-dimensional  translational  transverse  amplitude 
for  airfoil  1 

-  Non-dimensional  translational  transverse  amplitude 
for  airfoil  2 

-  Phase  angle  in  degree  between  the  chordwise  and 
transverse  translational  oscillation  with  the  latter 
as  reference  for  the  first  airfoil. 

-  Phase  angle  in  degree  between  the  chordwise  and 
transverse  translational  oscillation  with  the  latter 
as  reference  for  the  second  airfoil. 

Format  (5F10.  6,15)  -  1  data  card. 

-  Final  non-dimensional  time  to  terminate  unsteady  flow 
solution. 

-  Starting  time  step  for  non-osc.  motions  if  TADJ  “0.0 

-  No.  of  computational  steps  per  cycle  for 
harmonic  motion. 

-  Baseline  time  step  size  for  all  motions  if  TADJ  #  0. 0 

-  Tolerance  criterion  for  checking  the  convergence 
between  successive  iterations  of  (Uw)k  and  (Py)*- 

-  Factor  by  which  DTS  will  be  adjusted. 

-  Steady  lift  coefficient  for  the  single  airfoil  at  the 
specified  AOA 

-  Option  for  changing  the  unsteady  Kutta  condition  to 
satisfy  the  tangential  velocity  as  per  Geising's  case. 

-  0  equal  pressure  at  the  trailing  edge  panels. 

-  1  equal  tgt  velocities  at  the  trailing  edge  panels. 
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TK 

TKM1 


-  Time  step  t%. 

-  Time  step 

ALPHA(L)  -  Angle  of  attack  of  airfoil  L  at  time  t^. 

OHEGA(L)  -  Rotational  velocity  (positive  counter  clockwise)  at  time 
tfr  for  airfoil  L 

U(L)  -  Chordwise  translational  velocity  (positive  forward)  at 
time  tjf  for  airfoil  L 

V(L)  -  Transverse  translational  velocity  (positive  downward)  at 
time  tfc  for  airfoil  L 

NITR  -  Iteration  number. 

VXW(L)  -  Iterative  solution  of  ( U „)%  of  airfoil  L 

VYW(L)  -  Iterative  solution  of  of  airfoil  L 

UAKE(L)  -  Iterative  solution  of  shed  vorticity  panel  length 

of  airfoil  L 

TKETA(L)  -  Iterative  solution  of  shed  vorticity  panel  orientation 
of  airfoil  L 

GAMK(L)  -  Iterative  solution  of  the  strength  of  the  current 
vorticity  distribution  of  airfoil  L 

J  -  Panel  number. 

XI (J)  -  local  x-coordinate  of  the  midpoint  of  j panel. 

YI(J)  -  local  y-coordinate  of  the  midpoint  of  j*b  panel. 

X(J)  -  global  X-coordinate  of  the  midpoint  of  jth  panel. 

Y(J)  -  global  Y-coordinate  of  the  midpoint  of  jth  panel. 

Q( J)  ■  Strength  of  source  distribution  on  the  J **  panel. 

CP(J)  -  Pressure  coefficient  at  the  midpoint  of  the  j*b  panel. 

V(J)  -  Total  tangential  velocity  at  the  midpoint  of  the  jth  panel 
with  respect  to  the  moving  local  system. 

Figure  8.  List  of  Output  Variables. 
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VN(J)  -  Normal  velocity  at  the  midpoint  of  the  j*b  panel 
with  respect  to  the  moving  local  system. 

PHIK(J)  -  Potential  at  the  mid-point  of  the  j ^  panel  at 
the  current  time  step. 

PHI(J)  -  Potential  at  the  mid-point  of  the  jt^1  panel  at 
previous  time  step. 

INT6AMMA  -  Integral  of  the  disturbance  velocity  around  the  airfoil. 

CD(L)  -  Drag  coefficient  of  airfoil  L. 

CL(L)  -  Lift  coefficient  of  airfoil  L. 

CM(L)  -  Pitching  moment  coefficient  about  leading  edge 
of  airfoil  L. 

M  -  Trailing  wake  core  vortex  number. 

X1I(M)  -  X-coordinate  of  the  center  of  the  m core  vortex 
of  airfoil  1  with  respect  to  local  system. 

Y1I(M)  -  Y-coordinate  of  the  center  of  the  «c*  core  vortex 
of  airfoil  1  with  respect  to  local  system. 

X1(M)  -  X-coordinate  of  the  center  of  the  core  vortex 

of  airfoil  1  with  respect  to  global  system. 

Y1(M)  -  Y-coordinate  of  the  center  of  the  core  vortex 

of  airfoil  1  with  respect  to  global  system. 

X2I(M)  -  X-coordinate  of  the  center  of  the  a core  vortex 
of  airfoil  2  with  respect  to  local  system. 

Y2I(M)  -  Y-coordinate  of  the  center  of  the  core  vortex 

of  airfoil  2  with  respect  to  local  system. 

X2(M)  -  X-coordinate  of  the  center  of  the  core  vortex 

of  airfoil  2  with  respect  to  global  system. 

Y2(M)  -  Y-coordinate  of  the  center  of  the  core  vortex 

of  airfoil  2  with  respect  to  global  system. 

CIRC(M,L)-  Circulation  strength  of  the  a core  vortex  of 
airfoil  L 

Figure  8  (Cont'd) 
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V.  RESULTS  AND  DISCUSSION  OF  CASE-RUNS 


USPOTF2  is  primarily  written  as  a  follow-up  to  U2DIIF  for  the  single  airfoil.  All 
the  case-runs  with  the  exception  of  the  gust  case  will  be  presented.  The  approach  in  the 
gust  case  is  not  consistent  with  the  requirements  for  an  irrotational  flow  Held  and  will 
not  be  treated  in  this  report.  The  step  change  in  angle  of  attack  (AOA)  will  be  compared 
with  Giesing's  for  the  same  airfoil  set  at  the  final  AOA,  undergoing  an  impulsive  start 
from  rest.  As  there  exist  no  comparison  data  for  the  other  sub-cases,  the  results  will  thus 
not  be  as  extensive  as  the  step  input.  They  are  documented  here  for  the  sole  purpose 
of  illustrating  the  capability  of  the  code. 

A.  STEP  CHANGE  IN  ANGLE  OF  ATTACK 

1.  Case  Run  Definition 

Consider  a  first  case  of  two  airfoils  initially  at  zero  AOA  to  the  free  stream  Vm 
which  undergo  an  out-of-phase  step  change  in  AOA  (attv)  at  time  t„.  Consider  a  second 
case  of  two  airfoils  at  rest,  set  initially  at  an  AOA  of  oc„„  out-of-phase  with  one  another, 
given  an  impulsive  start  to  V„.  The  above  two  cases  are  equivalent  within  the  thin 
airfoil  approximation.  Case  1  is  computed  by  USPOTF2  and  case  2  is  obtained  by 
Giesing's  computational  analysis. 

2.  Differences  between  USPOTF2  and  Giesing's  code 

While  both  codes  use  an  extension  of  the  PANEL  method  to  solve  for  the  un¬ 
steady  potential  flow  solution  for  two  airfoils,  perfect  correlation  is  not  possible  for  se¬ 
veral  reasons.  The  reader  is  referred  to  References  (5)  and  [7]  for  a  detailed  description 
of  Giesing's  approach.  Some  basic  differences  are: 

•  USPOTF2  models  the  wake  vortex  sheet  by  a  wake  element  and  a  series  of  point 
vortices  shed  through  the  wrake  element.  This  follows  the  approach  of  Basu  and 
Hancock  with  the  necessary  assumptions  on  the  wake  element  characteristics. 
Giesing  treated  the  wake  vortex  sheet  to  comprise  of  a  distributed  line  vortex 
which  is  convected  at  the  local  fluid  velocity  at  the  vortex  location  assuming  that 
the  small  portion  of  the  wake  being  shed  does  not  contribute  to  the  convection  of 
the  wake. 

•  The  time  step  increment  in  L  SPOTF2  is  a  parameter  that  afreets  the  overrall  results 
of  the  unsteady  flow  in  that  by  having  small  time  step,  the  point  vortices  being  shed 
will  be  greater  in  number  but  w  eaker  in  strength,  while  for  bigger  time  step,  the 
point  vortices  become  fewer  but  have  effectively  stronger  strengths.  The  time  step 
increment  in  Giesing's  code  is  important  in  that  the  assumption  that  the  small 


portion  of  the  wake  being  shed  does  not  contribute  to  the  convection  of  the  wake 
is  only  exactly  true  in  the  limit  when  the  time  step  becomes  zero. 

•  Giesing  uses  the  Adam's  formula  of  varying  degree^  to  convect  the  wake  vortex 
sheet  while  USPOTF2  uses  only  a  predictor  algorithm. 

•  The  treatments  of  the  circulation  T(/)  are  different  for  both  codes.  USPOTF2  ex¬ 
tended  the  influence  coefficient  concept  to  the  unsteady  flow  regime  and  solved  for 
the  circulation  using  the  Kutta  condition  of  equations  3.33  and  3.34.  All  calcu¬ 
lations  w'ere  done  in  the  moving  frame  in  the  presence  of  the  unsteady  wake  for¬ 
mations.  Giesing  considered  the  circulation  to  be  a  combination  of  the 
quasi-steady  circulation  and  the  circulation  due  to  the  vortex  wakes.  From  Refer¬ 
ence  [  5] 

HO  -  r^/)  +  r^/)  (5.1) 

where  T £1)  is  the  circulation  required  to  satisfy  the  Kutta  condition  on  body  (/)  as 
it  moves  through  the  fluid  when  it  is  assumed  that  the  body  does  not  shed  any 
vorticity  and  r,(/)  is  that  circulation  required  to  satisfy  the  Kutta  condition  on  body 
(/)  as  it  travels  through  the  flow  field  generated  by  the  vortex  wakes  shed  by  the  two 
bodies  and  the  flow  field  generated  by  the  circulatory  flow  about  the  other  body. 

•  Giesing's  Kutta  condition  reqires  equal  tangential  velocities  at  the  upper  and  lower 
surface  panels  at  the  trailing  edges  while  the  Kutta  condition  of  USPOTF2  pre¬ 
scribes  equal  pressures.  In  order  to  have  a  meaningful  comparison,  LSPOTF2  in¬ 
cludes  an  option  for  equal  tangential  velocities. 


•  USPOTF2  requires  the  actual  computation  of  the  total  velocity  potential  from  a 
combination  of  an  onset  flow  potential,  disturbance  potential  and  an  assumption 
of  a  reference  potential.  Giesing  uses  the  technique  of  the  Douglas  Potential  Flow 
I  program  which  treats  the  potential  as  the  combination  of  the  quasi-steady  potential 

and  the  potential  due  to  the  vortex  wakes  where  these  two  terms  are  defined  as 
before.  The  total  velocity  potential,  in  his  case,  need  not  be  computed  but  only  the 
time  derivative  of  it  which  is  then  written  in  terms  of  previously  known  parameters. 

.  •  The  number  of  panels  and  its  distribution  on  the  airfoils  are  different  for  both 

I  codes.  This  will  be  accentuated  at  the  trailing  edges  and  will  lend  itself  to  differ¬ 

ences  in  the  trailing  edge  panel  distribution. 


3.  Results  and  Discussions 

Figure  9  shows  Giesing's  results  for  the  Von  Mises  8.4  per-cent  thick,  symmet¬ 
rical  airfoil  undergoing  an  impulsive  start.  This  can  be  compared  with  Figure  10  which 
shows  the  result  obtained  with  USPOTF2  for  the  same  time  step  with  NGIES  set  to 
onei’.  The  results  are  not  perfectly  correlated  but  the  quality  and  order  of  magnitude 
agreement  are  excellent.  Figure  11  gives  essentially  the  same  result  but  with  NGIES  set 


18  This  is  essentially  a  predictor-corrector  algorithm. 

19  This  Kutta  condition  results  in  equal  tangential  velocities  at  the  trailing  edge  panels  of  the 
airfoils. 
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to  zero20.  Surprisingly,  the  results  obtained  for  both  types  of  Kutta  condition  turn  out 
to  be  quite  similar.  This  is  seen  especially  in  the  pressure  coefficient  plots  when  the  two 
plots  are  put  together  as  seen  in  Figure  12.  The  pressure  coefficient  agrees  over  70  per¬ 
cent  of  chord  length  at  the  lower  surface  and  over  90  per-cent  of  chord  length  at  the 
upper  surface  with  the  greatest  discrepancy  at  the  trailing  edge.  Further  comparison  for 
smaller  angles  of  attack  would  be  of  interest  to  see  whether  this  similarity  is  generally 
true. 

Figure  13  compares  the  aerodynamic  characteristics  when  the  vertical  distance 
YSHIFT  is  varied.  Figure  14  gives  the  time  variation  for  a  larger  time  step  of  the 
normalised  lift,  moment  coefficient  and  drag  coefficient  for  the  particular  case  of 
YSHIFT  *  2.0. 


20  This  Kutta  condition  results  in  equal  pressure  coefficients  at  the  trailing  edge  panels  of  the 
airfoils. 


J:. ■-.1%.:.  ■ 
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Figure  9.  Giesing's  Calculated  Results:  Impulsive  Start  for  an  S.4  per-cent  thick 
Von  Mises  airfoil  set  at  a  *=  0.8  radians  for  YSHIFT  *  2.0  [reproduced 
with  permission  from  Ref.  5], 


i 
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(b)  Time  History  of  the  Lift  Coefficients. 


figure  9  (Cont'd) 


Figure  10.  USPOTF2  Results  obtained  with  USPOTF2  for  equal  velocities  at  the 
trailing  edge.:  Step  change  in  AOA  for  an  8.4  per-cent  thick  Von 
Mises  airfoils  placed  at  2  chord  length  vertical  distance  with  initial 
AOA  *  0.0  radians  and  final  AOA  *  0.8  radians  pivoting  at  the 
leading  edges. 
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(c)  Time  History  of  the  Lift  Coefficients. 


Figure  10  (Cont'd) 


(a)  Pressure  Distribution  at  tVJc 


Figure  1 1.  USPOTF2  Results  obtained  with  USPOTF2  for  equal  pressures  at  the 
trailing  edge.:  Step  change  in  AOA  for  an  8.4  per-cent  thick  Von 
Mises  airfoils  placed  at  2  chord  length  vertical  distance  with  initial 
AOA  *  0.0  radians  and  final  AOA  *  0.8  radians  pivoting  at  the 
leading  edges. 


(b)  Wake  pattern  at  tVoa/c  *  0.65 


Figure  11  (Cont'd) 
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(c)  Time  History  of  the  Lift  Coefficients. 

i  -  ' - 

Figure  1 1  (Cont'd) 
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8  :  Equal  Pressure  Coefficient  at  trailing  edge  panels. 
1  :  Equal  Tangential  Velocities  at  trailing  edge  panels 


0 


0.2 


OA 

NON— DIMENSIONAUSED  TIME 


(b)  Time  History  of  the  Lift  Coefficients. 


Figure  13  (Cont'd) 


NON— DIMENSI0NAL1SED  TIME 


(a)  Time  History  of  the  Lift  Coefficients. 

Figure  14.  Step  Change  in  Angle  Of  Attack:  Step  change  in  AOA  for  an  8.4 
per-cent  thick  Von  Mises  airfoils  placed  at  2  chord  length  vertical  dis¬ 
tance  with  initial  AOA  =  0.0  radians  and  final  AOA  *  0.8  radians 
pivoting  at  the  leading  edges. 


(b)  Time  History  of  the  Moment  Coefficients. 


Figure  14  (Cont'd) 


Figure  14  (Cont'd) 

B.  OTHER  SUB-CASES 

The  reader  is  referred  to  the  results  and  disussions  of  Reference  (  1  ]  for  the  com¬ 
parison  of  the  single  airfoil  case  with  existing  codes  in  the  phase  and  magnitude  re¬ 
lationship  of  the  aerodynamic  coefficients  to  the  forcing  function  as  well  as  the  wake 
convection.  The  present  results  show  that  the  same  trend  is  maintained  as  in  the  single 
airfoil  case  with  no  significant  phase  shift  but  with  a  general  magnitude  change  due  to 
an  effective  ground  effect  caused  by  the  presence  of  the  second  airfoil. 
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1.  Modified  Ramp  Change  in  Angie  of  Attack 

As  for  the  single  airfoil,  the  modified  ramp  is  defined  mathematically  as  follows: 

0  t  <  0 

a(f)  =  <5a(3  —  2r/r)/2/T2  0  £  /  <,  r  (5.2) 

da  t  >  t 

where  da  is  the  magnitude  of  the  AOA  change  and  r  is  the  rise  time  for  the  AOA  to 
reach  its  final  value.  Figure  15  treats  the  case  of  the  2  airfoils  undergoing  an  out-of¬ 
phase  modified  ramp  change  in  the  angle  of  attack  of  0.1  radians  with  a  time  constant 
of  1.5.  The  plots  show  the  effects  on  the  aerodynamic  coefficients  due  to  the  presence 
of  the  second  airfoil.  The  corresponding  aerodynamic  coefficients  of  the  single  airfoil 
are  plotted  for  comparison  purposes. 

2.  Translational  Harmonic  Motion 

The  code  is  capable  of  computing  the  unsteady  flow  solution  for  any  general 
translational  motion  described  by  a  chordwise  and  a  transverse  component  bearing  a 
given  phase  relationship  with  the  restriction  that  the  2  airfoils  move  only  in-phase  or 
out-of-phase.  The  translational  harmonic  motion  is  described  by 

hy{i)  =  dhy  Sin  (cot) 

hx(i)  =  dhx  Sin  (cor  +  /)  (5.3) 

where  to  is  the  oscillation  frequency,  /  is  the  phase  angle  between  the  chordwise  and 
transverse  oscillation  and  dh„  and  dht  are  the  magnitudes  of  chordw»se  and  transverse 
oscillations  respectively.  The  case-run  considered  in  this  section  relates  to  a  pure  heav¬ 
ing  or  plunging  motion.  A  NACA-0015  airfoil  is  chosen  for  the  case-run.  The  airfoils 
are  set  at  zero  radian  angles  of  attack  and  subsequently  given  an  out-of-  phase  plunging 
oscillation  at  an  amplitude  of  0.018  chord  length  at  a  non-dimensionalised  frequency  of 
1.  Figure  16  shows  the  aerodynamic  coefficients  for  the  two  airfoils  and  compares  them 
with  the  single  airfoil  case  where  applicable.  Note  that  the  plots  for  the  trailing  wake 
uses  different  scales  for  the  x  and  y  axes. 

3.  Rotational  Harmonic  Motion 

The  treatment  of  the  harmonic  pitching  motion  is  similar  to  the  modified  ramp 
case.  As  for  the  other  sub-cases,  the  airfoils  are  restricted  to  in-phase  and  out-of  phase 
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motion.  The  case  of  the  out-of-phase  motion  will  be  treated  here.  The  harmonic 
pitching  oscillation  is  described  by: 


a(t)  -  5a  Sin  (cor)  (5.4) 

where  5a  and  c o  are  the  amplitude  and  frequency  of  the  harmonic  oscillation  respec¬ 
tively.  Figure  17  shows  the  results  of  the  8.4%  thick  Von  Mises  symmetric  airfoil  os¬ 
cillating  at  an  amplitude  of  0.1  radian  at  a  reduced  frequency  of  c oc/Vx  =  20.0  about  the 
leading  edge.  Again  the  plots  are  given  together  with  the  single  airfoil  case  undergoing 
the  same  motion  for  comparison  purposes. 
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(a)  Time  History  of  the  Lift  Coefficients. 

Figure  IS.  Modified  Ramp  Change  in  Angle  Of  Attack:  Ramp  change  in  AOA 
for  an  8.4  per-cent  thick  Von  Mises  airfoils  placed  at  1  chord  length 
vertical  distance  with  initial  AOA  =  0.0  radians  and  final  AOA  =  0.1 
radians,  rise  time  of  1.5  pivoting  at  the  mid  chord. 


(b)  Time  History  of  the  Moment  Coefficients 


Figure  15  (Cont'd) 
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(c)  Time  History  of  the  Drag  Coefficients. 
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(d)  Time  History  of  the  Circulation 


Figure  15  (Cont'd) 
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(0  Wake  pattern  at  tVJc  ■  2.3 


Figure  15  (Cont'd) 
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(a)  Time  History  of  the  Lift  Coefficients. 

Figure  16.  Harmonic  Plunging  Motion:  Translational  harmonic  plunging  AOA 
for  \ACA-0015  airfoils  placed  at  1  chord  length  vertical  distance  set 
at  AOA  =  0.0  radians  with  plunging  amplitude  of  0.018  chord  length 
at  a  reduced  frequency  of  1. 
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(a)  Time  History  of  the  Lift  Coefficients. 

Figure  17.  Harmonic  Pitching  Motion:  Rotational  harmonic  pitching  AOA  for 
NACA-0015  airfoils  placed  at  1  chord  length  vertical  distance  set  at 
AOA  =  0.0  radians  with  pitching  amplitude  of  0.1  radian  at  a  reduced 
frequency  of  4  pivoting  about  the  leading  edges. 
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(c)  time  History  of  the  Drag  Coefficients 


Figure  17  (Cont'd) 
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(d)  Time  History  of  the  Circulation. 
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Figure  17  (Cont'd) 


(e)  Pressure  Coefficient  at  tvm  ■  3,0. 


Figure  17  (Cont'd) 
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VI.  CONCLUSION 


A.  GENERAL  COMMENTS 

LSPOTF2  has  been  developed  as  an  intermediate  stage  to  getting  a  solution  for  a 
full  cascade  undergoing  unsteady  motion.  In  itself,  it  can  be  used  to  simulate  unsteady 
wing-tail,  wing-canard,  wing-flap,  wing-aileron,  horizontal  tail-elevator,  vertical 
tail-rudder  and  a  whole  host  of  moving-stationary  airfoil  interactions.  These  simu¬ 
lations  require  a  little  addition  to  the  main  program  and  would  appear  as  modules 
modelling  the  variation  of  the  global  angles  of  attack  and  relative  displacements  as  per 
all  the  sub-cases  done.  The  subroutines  would  not  be  affected  by  the  severity  of  the  test 
case  except  for  possibly  subroutine  NEWPOS. 

Validation  of  LSPOTF2  has  been  done  against  Giesing's  code  for  the  particular  case 
of  a  step  change  in  angle  of  attack.  This  was  shown  to  have  good  correlation.  However, 
it  can  not  be  said  that  the  agreement  will  hold  for  other  angles  of  attack  and  more  sub¬ 
case  runs  for  smaller  angles  of  attack  should  be  made.  In  the  same  manner,  the  success 
of  the  step  change  in  AOA  in  no  way  validates  the  other  sub-cases  which,  for  the  time 
being,  remain  unproven. 

B.  ENHANCING  USPOTF2  PROGRAM  CAPABILITY 

As  noted  above,  LSPOTF2  needs  to  be  tried  more  extensively  either  with  the  exist¬ 
ing  sub-cases  or  with  new  sub-cases  against  existing  numerical  or  experimental  results. 

The  software  in  itself  has  an  implicit  weakness  in  the  numerical  computation  of  the 
velocity  potential.  The  velocity  potential  at  the  leading  edge  is  obtained  by  integrating 
the  velocity  field  from  the  leading  edge  to  a  point  100  chord  lengths  upstream  of  the 
leading  edge.  Two  assumptions  are  made:  First,  the  disturbance  velocity  at  the  point 
100  chord  lengths  of  the  leading  edge  approaches  zero.  Second,  the  contribution  to  the 
velocity  potential  due  to  the  integration  from  the  point  100  chord  lengths  upstream  of 
the  leading  edge  to  upstream  infinity  must  not  change  in  time.  The  coefficients  of  pres¬ 
sure  will  be  accurate  only,  if  those  two  assumptions  hold  (at  least  in  an  approximate 
sense).  This  however  does  not  give  a  very  exact'  solution  as  by  increasing  the  chord 
length  by  900  per-cent  to  1000  chord  lengths,  there  is  a  change  in  the  pressure  coefficient 
at  the  trailing  edge  of  about  8  per-cent  for  the  first  time  step.  Krainer  has  implemented 
an  analytical  solution  to  the  potential  problem  for  the  single  airfoil.  A  similar  procedure 
to  upgrade  USPOTF2  should  be  considered. 
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Other  improvements  to  the  code  would  be  the  reduction  of  redundant  computations. 
An  obvious  example  would  be  the  subroutines  VELD1S  and  PRESS  which  essentially 
do  the  same  work  for  the  steady  and  unsteady  flow  respectively.  Again  Krainer  [Ref. 
6]  has  improved  the  original  code  for  the  single  airfoil  (U2DI1F)  and  though  some  of 
his  improvements  were  implemented  in  USPOTF2,  there  still  remains  a  task  to  do  the 
full  job  completely. 

Finally,  the  original  primary  objective  of  extending  the  code  to  solve  for  the  un¬ 
steady  flow'  solution  for  3,4  airfoils  and  leading  to  a  full  cascade  still  remains  a  big  task, 
not  just  for  the  programmer  but  also  for  the  computer  system  in  terms  of  computational 
storage  and  time  requirements. 
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APPENDIX  A.  USPOTF2  SOURCE  LISTINGS 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PROGRAM  USP0TF2 
VERSION  1 
SEPTEMBER  88 

UNSTEADY  MOTION  FOR  TWO  AIRFOILS  IN  POTENTIAL 
INCOMPRESSIBLE  FLOW 

USING  PANEL  METHODS  BASED  ON  THE  HESS  &  SMITH 
AND  WAKE  MODEL  BASED  ON  BASU  &  HANCOCK 

icccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y( 202) , 

+  C0STHE( 201) ,SINTHE( 201) ,SS(2) ,NP1 ,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRF0>XI(202) ,YI(202) , 

+  C0STHL( 20 1 ) , S INTHL( 201) 

COMMON  /WAK/  VYW(2) ,VXW( 2) ,WAKE( 2) ,DT 
COMMON  /WAK2/  VYVK(2) ,VXWK(2) 

COMMON  /SING/  Q( 200) ,GAMMA(2) ,QK( 200) ,GAMK(2) 

COMMON  /CORV/  CV(2,200)  ,XC(2,200)  ,YC(2,200)  ,M,TD,CWX(2,200)  , 

+  CWY(2 ,200)  ,XCI(2,200)  ,  YCI(2 ,200) 

COMMON  /POT/  PHI(200) ,PHIK(200) 

COMMON  /GUST/  UG( 200 ) , VG( 200 ) , XGF , UGUST , VGUST 
COMMON  /EXTV/  UE(200) ,VN(200) 

COMMON  /PARD/  NACAD( 2) ,TAUD(2) ,EPSMAD(2) ,PTMAXD(2) 

COMMON  /COF/  A( 201, 211) ,KEQNS 

COMMON  /GEOM/  SINALF(2) ,COSALF( 2) ,OMEGA( 2) ,UX(2) ,UY(2) , PIVOT, 

+  XPRM.YPRM 

COMMON  /CPD/  CP( 200) ,SCL,T,SCM,SGAM 
COMMON  /GIES/  NGIES 

DIMENSION  XXC( 2,200) , YYC( 2,200) ,T0L1( 2) ,TOL2( 2) ,THENP1( 2) , 
+SINANG(2) ,COSANG(2) ,C0SDA(2) ,SINDA(2) ,DHX(2) ,DHY(2) ,ALP(2) , 

+ALPHA( 2 ) , ALP I ( 2 ) , ANGLE ( 2 ) , DELHX( 2 ) , DELHYC 2 ) , PHASE( 2 ) , PHA( 2 ) , 

+HX(2) , HXO( 2 ) ,HY(2) ,HY0(2) 

INITIALISATION 

PI  =  3. 1415926585 

M  =  0 

XPRM  =0.0 

YPRM  =  0.0 

INPUT  FROM  FILE  CODE  5  AND  SET  UP  LOCAL  PANEL  NODES  AND  SLOPES 

WRITE  (6,1003) 

1003  FORMAT  (/////,'  DATA  READ  FROM  FILE  CODE  1 ',//) 

CALL  INDATA 
CALL  SETUP 

READ  (1,502)  ALPI(1),ALPI(2),DALP,TC0N, FREQ, PIVOT, IPHASE 
WRITE  (6,502)  ALPI( 1) ,ALPI( 2) ,DALP,TCON,FREQ, PIVOT, IPHASE 
501  FORMAT  (8F10. 6) 


onoocoooooo 


o  n  n  n  o 


■EJUW'JWS? 


502  FORMAT  (6F10.6.I3) 

READ  (1,501)  UGUST ,VGUST,DELHX( 1) ,DELHX( 2) ,DELHY( 1) ,DELHY(2) , 
+  PHASE ( 1 ) , PHASE ( 2 ) 

WRITE  (6,501)  UGUST, VGUST,DELHX( 1) ,DELHX( 2) ,DELHY( 1) ,DELHY(2) , 
+  PHASE ( 1) ,PHASE( 2) 


READ  (1,503)  TF , DTS ,TOL ,TAD J , SCL , SCM , SGAM , NGIES 
WRITE  (6,503)  TF, DTS, TOL, TAD J, SCL, SCM, SGAM, NGIES 
503  FORMAT  (7F10.6.I3) 

IF  (IFLAG  .EQ.  0)  WRITE  (6,1005) 

1005  FORMAT  (///,'  COORDINATES  OF  AIRFOIL  NODES', 

+  // ,3X,  X/C* ,6X, '  Y/C',/) 

IF  (IFLAG  .EQ.  0)  WRITE  (6,1010)  (XI(I) ,YI(I) ,I=l,2*(NODTOT+l)) 
1010  FORMAT  (F10.  6,F10.  6) 

WRITE  (6,1000)  NAIRFO 

1000  FORMAT(//,‘ TOTAL  NO  OF  AIRFOILS  -  ’,14,/) 


STEADY  FLOW  CALCULATION  AT  ALPI(I) 


INITIAL  DEFINITION 


NP1  =  NODTOT  +  1 

NP2  *  2  *  NP1  +  1 

NP3  ■  NAIRF0*N0DT0T+1 

NP4  =  NP3+1 

NP5  *  NP4+1 

DO  101  I  =  1,2 
ALPHA(I)  *  ALPI(I) 

IF  ( ALPHA ( I )  .GT.  90.)  GO  TO  200 
COSALF(I)  =  COS(ALPHA( I)*PI/180.  ) 

101  SINALF(I)  *  SIN(ALPHA(I)*PI/180.  ) 

CALL  NEWPOS(O) 

DO  1100  L  «  1, NAIRFO 
1100  WRITE  (6,1020)  L,SS(L) 

WRITE  (6,1040)  XSHIFT.YSHIFT 

1020  FORMAT(//, '  AIRFOIL( '  12, ' )  PERIMETER  LENGTH  ■  '.F10.6,/) 

1040  FORMAT( // , '  XSHIFT  *  ,F5. 1,'  YSHIFT  ®  ' ,F5. 1,/) 

DO  102  1-1,2 

102  WRITE  (6,1030)  I,ALPHA(I) 

1030  FORMAT  (//,'  STEADY  FLOW  SOLUTION  AT  ALPHA( ' ,12, ' )  «  ’,F10.6,/) 
IF  (SCL. NE. 0)  WRITE  (8,1035) 

1035  FORMAT( 3X , ’ TIME ’ , 4X , ’ CL( 1 ) /SCL ' , IX , ' CL( 2 ) /SCL ' , IX , ’ CM( 1 ) /SCM ' , 
+1X,'CM(2)/SCM' ,3X, ' CD( 1) ' ,4X,,CD(2)' ,4X, 
+'GAMK(1)/SG,,1X,'GAMK(2)/SG' ,//) 

IF  (SCL.  EQ. 0)  WRITE  (8,1036) 

1036  F0RMAT(3X  'TIME' ,6X, *CL( 1) ' ,5X, ’CL(2) ' ,5X, 'CM( 1) ' , 

+5X, ’CM(2)' ,5X, 'CD(l)' ,4X, 'CD(2)'//) 

CALL  I NFL  (0) 

CALL  COEF  (0) 

CALL  GAUSS(3,0,0) 

CALL  KUTTA( NITR , PVTAG) 

CALL  VELD IS 
SIN1  =  SINALF(l) 

SIN2  -  SINALF(2) 

C0S1  =  COSALF(l) 

C0S2  =  C0SALF(2) 
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CALL  FANDM(SIN1,SIN2,C0S1,C0S2) 

INITIALISATION  FOR  UNSTEADY  FLOW  CALCULATION  TO  BEGIN 

DA  =0.0 
XGF  =  0.  0 
DO  100  L  =  1.NAIRF0 
HX(L)  =0.0 
HY(L)  =0.0 
HXO(L)  =0.0 
HYO(L)  =0.0 

ANGLE(L)  =  ALPI(L)*PI/180.  +  ATAN( VGUST/C 1.  +UGUST) ) 
OMEGA(L)  =0.0 
ALP(L)  =  ALPI(L) 

DHX(L)  =0.0 

DHY(L)  =0.0 

COSANG(L)  =  COS(ANGLE(L)) 

SINANG(L)  =  SIN(ANGLECL)) 

UX(L)  =  0.  0 

UY(L)  =  0.  0 

COSDA(L)  =1.0 
SINDA(L)  =0.0 
KIG  =  (L  -1)*N0DT0T 

VXW(L)  =  COSALF(L) 

VYV(L)  =  SINALF(L) 

GAMK(L)  =  GAMMA(L) 

PHA(L)  =  PHASE(L)*PI/180. 

DO  100  IG  =  1 ,  NODTOT 
UG( IG+KIG)  =0.0 
100  VG( IG+KIG)  =0.0 
T  =0.0 

TOLD  =0.0 

RIGID  BODY  MOTIONS  OF  AIRFOIL 

IF  (FREQ  .NE.  0.0)  GO  TO  1 
IF  (DALP  .EQ.  0.0)  GO  TO  2 
IF  (TCON  .NE.  0.0)  GO  TO  3 
IF  ( I PHASE  .EQ.  0)  GO  TO  4 
ALPHA( 1)=  ALPI(l)  +  DALP 
ALPHA( 2)=  ALPI(2)  +  DALP 
GO  TO  5 

4  ALPHA( 1)=  ALPI(l)  +  DALP 
ALPHA(2)=  ALPI(2)  -  DALP 

5  DO  6  I  =  1,2 

COSALF(I)  =  C0S(ALPHA(I)*PI/180. ) 

SINALF(I)  =  SIN(ALPHA(I)*PI/180.  ) 

6  CONTINUE 
CALL  NEWPOS(O) 

3  DT  =  DTS 

TD  =  DTS 

GO  TO  60 

2  IF  ((UGUST  .EQ.  0.0)  .AND.  (VGUST  .EQ.  0.0))  GO  TO  200 
DT  =  DTS 

TD  =  DTS 

GO  TO  60 
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1  OT  »  2. 0*PI/(FREQ*DTS) 

ID  *  DT 

60  T  *  DT 

WRITE  (6,1051) 

105 1  FORMAT  (/////»'  *************************************** '  / 

+  '  ***  begin’ UNSTEADY  FLOW  SOLUTION  ****',/, 

4-  *  it**********'*1***'***"***'****************** '  ) 

TRY  =0 

40  M  ■  M  +  1 

IF  (T  .GT.  TF)  GO  TO  200 

STORE  CORE  VORTEX  COORDINATES  FOR  TIME  STEP  ADJUSTMENTS 

IF  (M  .EQ.  1)  GO  TO  50 
DO  51  L  *  l.NAIRFO 

DO  51  I  *  l.M-l 

XXC(L.I)  «  XCI(L,I) 

YYC(L.I)  =  YCI(L, I) 

IF  (FREQ  .NE.  0.0)  GO  TO  11 

IF  (DALP  .EQ.  0.  0)  GO  TO  22 

IF  (TCON  -NE.  0.  0)  GO  TO  32 

STEP  CHANGE  IN  AOA 

IF  (TADJ  -NE.  0.0)  GO  TO  70 

IF  INCREMENTAL  PROGRESSIVE  TIME  STEP  IS  REQUIRED 
USE  . . .  TD  *  FLOAT( M+ 1 )*DTS  . . .  OTHERWISE  CONSTANT  TIME  STEP 
TD  *  DTS 

GO  TO  70 

MODIFIED  RAMP  CHANGE  IN  AOA 

33  IF  (T  .GT.  TCON)  GO  TO  34 

DAL  *  DALP  *  (3.  -2.  *T/TCON)*(T/TCON)**2 

OMEGA(l)  =  -  (DALP*PI/180. )  *  (6.  *T/(TCON*TCON))  *  (1. -T/TCON) 

IF  ( I PHASE  .EQ.  0)  GO  TO  41 

ALPHA(1)=  ALPI(l)  +  DAL 

ALPHA(2)=  ALPK2)  +  DAL 

OMEGA(2)  =  OMEGA(l) 

GO  TO  42 

41  ALPHA(1)=  ALPI(l)  +  DAL 

ALPHA(2)«  ALPI(2)  -  DAL 
OMEGA( 2 )  =  -  OMEGA(l) 

42  DO  43  I  *  1,2 

COSALF(I)  «  C0S(ALPHA(I)*PI/180.  ) 

SINALF(I)  «  SIN(ALPHA(I)*PI/180.  ) 

DA  «  ALPHA(I)  -  ALP(I) 

COSDA(I)  «  C0S(DA*PI/180. ) 

SINDA(I)  «  SIN(DA*PI/180. ) 

DHX(I)  «  PIVOT  *  (1. -COSDA(I)) 

DHY(I)  ■  -  PIVOT  *  SINDA(I) 

UY(I)  »  PIVOT  *  OMEGA(I) 

43  CONTINUE 
MTCON  =  M 
CALL  NEWPOS(O) 

GO  TO  70 
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34  DAL  =0.0 

IF  ( I PHASE  .EQ.  0)  GO  TO  45 
ALPHA( 1)=  ALPI(l)  +  DALP 
ALPHA(2)=  ALPI(2)  +  DALP 
GO  TO  46 

45  ALPHA(l)  =  ALPI(l)  +  DALP 

ALPHA(2)  =  ALPI(2)  -  DALP 

46  DO  44  I  =  1,2 

COSALF(I)  =  COS(ALPHA(I)*PI/180. ) 

SINALF(I)  =  SIN(ALPHA( I)*PI/180. ) 

DA  =0.0 

COSDA(I)  =1.0 
SINDA(I)  =0.0 
OMEGA(I)  =0.0 
DHX(I)  =0.0 
DHY(I)  =0.0 

UX(I)  =0.0 

UY(  I )  =  0.  0 

44  CONTINUE 

CALL  NEWPOS(O) 

IF  (TADJ  .NE.  0.0)  GO  TO  70 
TD  =  FLOAT( M+ 1 - MTC ON ) *DTS 

GO  TO  70 

SHARP  EDGE  GUST  (UGUST  AND/OR  VGUST)  NOT  PROVEN  DUE  TO 
INCONSISTENT  ASSUMPTIONS  REQUIRING  ROTATIONAL  FLOW 


22  XGF  =  T 

DO  113  L  =  l.NAIRFO 

LIG  =  (L-1)*N0DT0T 

KIG  =  (L-1)*NP1 

DO  110  IG  =  1 , NODTOT 
UG( IG+LIG)  =0.0 
VG(  IG+LIG)  =0.0 
XG  =  XCIG+KIG) 

XGP1  =  X( IG+KIG+1) 

IF  (IG  .LT.  NLOWER+1)  GO  TO  120 
IF  (XGF  .LE.  XG)  GO  TO  110 
IF  (XGF  .GE.  XGI1)  GO  TO  111 
FAC  =  (XOF  -  XG)/(XGP1  -  XG) 
UG(  IG+LIG)  =  UGUST*FAC 
VG( IG+LIG)  =  VGUST+FAC 
GO  TO  110 

111  UG( IG+LIG)  =  UGUST 
VG( IG+LIG)  =  VGUST 
GO  TO  110 

120  IF  (XGF  .LE.  XGP1)  GO  TO  110 
IF  (XGF  .GE.  XG)  GO  TO  121 

FAC  =  (XGF  -  X0P1)/(XG  -  XGP1) 
UG( IG+LIG)  =  UGUST*FAC 
VG( IG+LIG)  =  VGUST+FAC 
GO  TO  110 

121  UG( IG+LIG)  =  UGUST 
VG( IG+LIG)  =  VGUST 

110  CONTINUE 
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113  CONTINUE 

IF  (XGF  .LE.  COSALF(L))  MGUST  =  M 
IF  (TADJ  .NE.  O.O)  GO  TO  70 

IF  (XGF  .GT.  COSALF(L) )  TD  *  FLOAT(M+l-MGUST)*DTS 
GO  TO  70 
C 

C  TRANSLATION  HARMONIC  OSCILLATION 

C 

11  IF  (DALP  .NE.  0.0)  GO  TO  12 
DO  131  I  =  l.NAIRFO 

HX( I)  ■  DELHX(I)  *  SIN(FREQ*T  +  PHA(I)) 

HY(I)  =  DELHY(I)  *  SIN(FREQ*T) 

DHX(I)  «  HX( I)  -  HXO(I) 

DHY(I)  *  HY( I)  -  HYO(I) 

UX(I)  *  DELHX(I)*FREQ*COS(FREQ*T+PHA(I)) 

UY(I)  «  DELHY( I)*FREQ*COS(FREQ*T) 

131  CONTINUE 

XPRM  =  DHX( 2)  -  DHX(l) 

YPRM  *  DHY(2)  -  DHY(l) 

CALL  NEWPOS(O) 

GO  TO  70 
C 

C  ROTATIONAL  HARMONIC  OSCILLATION 

C 

12  DAL  «  DALP*SIN(FREQ*T) 

OMEGA(l)  =  -  (DALP*PI/180.  )  *  FREQ  *  COS(FREQ*T) 

IF  ( IPHASE  .EQ.  0)  GO  TO  141 
ALPHA(1)=  ALPI(l)  +  DAL 
ALPHA(2)=  ALPI(2)  +  DAL 
OMEGA( 2)  *  OMEGA(l) 

GO  TO  142 

141  ALPHA(1)=  ALPI(l)  +  DAL 

ALPHA(2)=  ALPI(2)  -  DAL 
OMEGA( 2 )  *  -  OMEGA ( 1 ) 

142  DO  143  I  *  1.2 

COSALF ( I )  *  COS(ALPHA(I)*PI/180.  ) 

SINALF(I)  *  SIN(ALPHA(I)*PI/180.  ) 

DA  *  ALPHA(I)  -  ALP(I) 

COSDA(I)  =  COS(DA*PI/180. ) 

SINDA(I)  «  SIN(DA*PI/180. ) 

DHX(I)  =  PIVOT  *  (1. -COSDA(I)) 

DHY(I)  *  -  PIVOT  *  SINDA(I) 

UY(I)  =  PIVOT  *  OMEGA(I) 

143  CONTINUE 
CALL  NEWPOS(O) 

C 

C  TRANSFORM  CORE  VORTEX  COORDINATES  V.  R.  T.  NEW  AIRFOIL  POSITION 

u 

70  IF  (M  .EQ.  1)  GO  TO  80 
DO  85  L  «  l.NAIRFO 

DO  90  I  *  l.M-l 

XCI(L.I)  *  XXC(L.I)  +  CWX(L.I)  *  DT 
YCI(L.I)  *  YYC(L.I)  +  CWY(L.I)  *  DT 
XCO  «  XCI(L.I) 

YCO  *  YCI(L.I) 

XCI(L.I)  *  XCO*COSDA(L)  -  YCO*SINDA(L)  +  DHX(L) 


85  CONTINUE  =  XC0*SINDA(L)  +  YCO*COSDA(L)  +  DHY(L) 

CALL  NEWPOS(3) 

80  CONTINUE 

WRITE  (6,1001)  T,DT 

1001  FORMAT  (/////,'  TIME  STEP  TK  =  ’ ,F10.  6 , 10X, 'TK  -  TKM1  =  '  F10  6  /) 
VKITZ  (6,1004)  ALPHA(  1)  , ALPHA(2)  !oMEGA( 1)  .OMEGAU^ D  ,UX(2) , ' ' ^ 

1004  FORMAT  (/,'  ALPHA(l)  =  '  F10.  6,5X,'  ALPHA(2)  a  ' ,F10. 6,/ 

+  OMEGA(l)  =  ,F10.6,5X,*  OMEGA(2)  «',F10.6, 

+/.’  0(1)  -  ’,F10.6,5X,’  U(2) 

+F10. 6,/,  V(l)  ■  1 ,F10. 6,5X,  V(2)  ■  1 ,  F10  6  III 

+  IX,  NITR  VXW(l)  VYV(l)  WAKE(l)  THETA(l)  GAMK(l) 

„  +  VXW(2)  VYV(2)  WAKE(2)  THETA(2)  GAMK(2)  ’,/) 

L» 

l  CALCULATE  THE  TRAILING  EDGE  WAKE  ELEMENT 

NUM  =0 

NITR  =  0 

10  DO  15  L  =  l.NAIRFO 

KI  =  (L-1)*(N0DT0T+1) 

WAKE(L)  *  SQRT(VYW(L)*VYW(L)+VXW(L)*VXW(L))*DT 
THENPl(L)  =  ATAN2(VYW(L),VXW(L)) 

COSTHL( NP1+KI )  »  COS(THENPKL)) 

15  SINTHL( NP1+KI)  =  SIN(THENP1(L) ) 

WRITE  (6,1002)  NITR,VXW( 1) ,VYW( 1) ,WAKE( 1) ,THENP1( 1) ,GAMK( 1) 

+VXW( 2) ,VYW( 2) ,WAKE( 2) ,THENP1( 2) ,GAMK( 2) 

1002  FORMAT  (15 ,4F10. 6,E14. 6,4F10. 6,E14. 6) 

XI(NP2)  *  XI(NPl)  +  WAKE(1)*C0STHL(NP1) 

YI(NP2)  =  YI(NPl)  +  WAKE(1)*SINTHL(NP1) 

XI(NP2+1)  *  XI(NPl)  +  WAKE ( 2)*C0STHL( 2*NP1 ) 

YI(NP2+1)  *  YI(NPl)  +  WAKE( 2)*SINTHL( 2*NP1) 

CALL  NEWPOS(l) 

CALL  INFL  (NITR) 

CALL  COEF  (NITR) 

CALL  GAUSS(3,M,NITR) 

WRITE  (6,*)  A(I,NP)',(A(I,101),I-1,100) 

WRITE  (6,*)  A(I ,NTOT) ' , ( A( I , 102) , I«1 , 100) 

WRITE  (6,*)' BEFORE  KUTTA' 

CALL  KUTTA( NITR , PVTAG ) 

IF  (PVTAG  .LT.  0.01)  GOTO  13 
NUM  *  NUM  +  1 
WRITE  (6,*) 'NUM  hum 
IF  (NUM  .  GT.  1)  GOTO  13 
DO  7  I  a  i,2 
VXW(I)  »  1 
7  VYW(I)  a  o 
GOTO  10 

13  CALL  TEWAK 

DO  24  L  «  l.NAIRFO 

T0L1(L)  »  ABS(VYW(L)  -  VYWK(L))/VYWK(L) 

T0L2(L)  =  ABS(VXW(L)  -  VXWK( L) ) / VXWK( L) 

VYW(L)  »  VYWK(L) 

24  VXW(L)  -  VXWK(L) 

IF  ((TOLl(l)  .LT.  TOL)  .AND.  (T0L2(1)  .  LT.  TOL) 

+.AND.  (T0L1(2)  .LT.  TOL)  .AND.  (TOL2(2)  .  LT.  TOL))  GO  TO  20 


98 


IF  (NITR  .  GT.  25)  THEN 
TOL  *  TOL* 10 

WRITE  (6,1023)  TOL 
ENDIF 

IF  (NITR  .GT.  50)  STOP 
NITR  *  NITR  +  1 

GO  TO  10 

20  WRITE  (6,1011)  NITR 

1011  FORMAT  (//,'  CONVERGED  SOLUTION  OBTAINED  AFTER  NITR  «  ’,13) 

1023  FORMAT  (//,'  *****  TOLERANCE  CRITERIA  CHANGED  :  TOL  «  ' ,F10. 6) 
CALL  PRESS 

IF  ((UGUST  .EQ.  0.0)  .AND.  (VGUST  .  EQ.  0.0))  GO  TO  300 
SIN1  «  SINANG(l) 

SIN2  =  SINANG(2) 

C0S1  «  COSANG ( 1) 

C0S2  *  COSANGf  2) 

CALL  FANDM( SIN1 , SIN2 , COS 1 , C0S2 ) 

GO  TO  400 

300  SIN1  »  SINALF(l) 

SIN2  ■  SINALF(2) 

C0S1  ■  COSALF(l) 

COS 2  =  COSALF(2) 

CALL  FANDM( SIN1 , SIN2 , C0S1 , COS2) 

400  CONTINUE 
C 

C  ADJUST  TIME  STEP  (TADJ  .  NE.  0.0)  IF  NECESSARY 
C 

IF  (TADJ  .EQ.  0.0)  GO  TO  95 
WRITE  (5,2001) 

2001  FORMAT  (//,'  DO  YOU  WANT  TO  ADJUST  TIME  STEP  ?  0  -  NO,  1  -  YES’) 
READ  (5,*)  IDT 
IF  (IDT  .EQ.  0)  GO  TO  95 
DT  «  TADJ  *  DT 

T  «  TOLD  +  DT 

WRITE  (6,1006) 

1006  FORMAT  (//,'  BACK -TRACK  COMPUTATION  AND  ADJUST  TIME-STEP',//) 

C 

C  WAKE  ELEMENT  LEAVES  TRAILING  EDGE  AS  A  CORE-VORTEX 

C 

95  DO  96  L  *  l.NAIRFO 

CV(L,M)  =  SS( L)*( GAMMA( L) -GAMK( L) ) 

CWX(L,M)  »  VXW(L) 

96  CWY(L.M)  «  VYW(L) 

XCI( 1,M)  «  XI(NPl)  +  0.5*WAKE(1)*C0STHL(NP1) 

YCI(l.M)  «  YI(NPl)  +  0.5*VAKE(1)*SINTHL(NP1) 

XCI(2,M)  ■  XI(NPl)  +  0. 5*WAKE(2)*C0STHL(NP1*2) 

YCI(2,M)  «  YI(NPl)  +  0.  5*WAKE(2)*SINTHL(NP1*2) 

CALL  NEWPOS(2) 

WRITE  (6,1052) 

1052  FORMAT  (//,'  TRAILING  VORTICES  DATA',//, 

+  4X,'M’  .AX.'Xl^M)'  .SX.’YIKM)’  ,6X. 'Xl(M)'  ,6X,'Y1(M)’ , 

+  5X, ’CIRCr ,6X, ’X2I(M) ' ,5X  ’ Y2I(M) f , 

+  5X,,X2(M)’  ,6X,'Y2(M)\6X,fCIRC2',/) 

DO  900  I  »  1,M 

900  WRITE  (6,1050)  I,XCI(1,I),YCI(1,I),XC(1,I),YC(1,I),CV(1,I), 

+  XCI(2, I) ,YCI(2, I) ,XC(2,I) ,YC(2,I) ,CV( 2,1) 


nnnnnnnonnnnonnnnon  on o 


1050  F0RMAT( 15 , 10F11. 6) 
CALL  C0RV0R 


RE -INITIALISE  PARAMETERS  FOR  NEXT  TIME  STEP  CALCULATION 

DO  30  L»  l.NAIRFO 
LI  =  (L-l)*NODTOT 
HXO(L)  «  HX(L) 

HYO(L)  =  HY(L) 

GAMMA(L)=  GAMK(L) 

ALP(L)  «  ALPHA(L) 

DO  30  I  «  l.NODTOT 
PHKI+LI)  «  PHIK(I+LI) 

30  CONTINUE 

TOLD  =  T 
DT  a  TD 
T  «  T  +  TD 

GO  TO  40 

200  WRITE  (6,1024)  TOL 

1024  FORMAT  (IX,/'*****  TOLERANCE  CRITERION  USED:  TOL  ' ,F10. 6) 

STOP 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  INDATA 


SET  PARAMETERS  OF  BODY  SHAPE 
FLOW  SITUATION,  AND  NODE  DISTRIBUTION 

USER  MUST  INPUT 

NLOWER  a  NUMBER  OF  NODES  ON  LOWER  SURFACE 
NUPPER  »  NUMBER  OF  NODES  ON  UPPER  SURFACE 
PLUS  DATA  ON  BODY  AND  SUBROUTINE  BODY 

SCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


NACAD(I)  .... 
TAUD(I)  .... 
EPSMAD(I)  ... 
PTMAXD(I)  ... 


NACA  NUMBERS  FOR  THE  TWO  AIRFOILS 

MAX  THICKNESS  FOR  THE  WO  AIRFOILS 

MAX  CAMBER  FOR  THE  WO  AIRFOILS 

CHORDWISE  POSN  FOR  MAX  CAMBER  FOR  THE  WO  AIRFOILS 


SUBROUTINE  INDATA 
DIMENSION  TITLE(20) 

COMMON  /BOD/  I FLAG, NLOWER, NUPPER, NODTOT,X( 202 ) ,Y( 202) , 

+  COSTHE(201) ,SINTHE(201) ,SS(2) ,NP1,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRFO,XI( 202) , YI( 202) , 

+  COSTHL(201) ,SINTHL(201) 

COMMON  /PARD/  NACAD(2) ,TAUD(2) ,EPSMAD(2) ,PTMAXD(2) 

READ  (1,501)  ITITLE 
WRITE  (6,501)  ITITLE 
DO  10  I  «  1, ITITLE 
READ  (1,502)  TITLE 
10  WRITE  (6,503)  TITLE 

501  FORMAT(3I5) 

502  FORMAT(20A4) 

503  FORMAT( IX , 20A4 ) 


100 


ooooooooooo 


onooonoooooooonoo  noon  non 


504  F0RMAT(I5,2F10.6) 

READ  (1,501)  IFLAG , NLOWER , NUPPER 
WRITE  (6,501)  IFLAG, NLOWER, NUPPER 
READ  (1,504)  NAIRFO.XSHIFT.YSHIFT 
WRITE  (6,504)  NAIRFO.XHIFT.YSHIFT 
IF  (IFLAG  . NE.  0)  RETURN 

COMPUTATION  FOR  NACA  4-DIGITS  SERIES 

READ  (1,501)  NACAD(l) 

WRITE  (6,501)  NACAD(l) 

READ  (1,501)  NACAD(2) 

WRITE  (6,501)  NACAD(2) 

DO  100  1=1 ,NAIRF0 

IEPS  «  NACAD( I)/1000 

IPTMAX  «  NACAD( I ) / 100  -  10*IEPS 

ITAU  «  NACAD(I)  -  1000*IEPS  -  100*IPTMAX 

EPSMAD(I)  «  IEPS*0. 01 

PTMAXD(I)  «  IPTMAX*0. 1 

TAUD(I)  «  ITAU*0.  01 

IF  (IEPS  .LT.  10)  GOTO  100 

COMPUTATION  FOR  NACA  5 -DIGITS  SERIES  NOTING  THAT  TAUD(I)  IS 
COMPUTED  AS  PER  4-DIGITS  SERIES. 

PTMAXD(I)  -0.2025 
EPSMAD(I)  «  2. 6595*PTMAXD(I)**3 
100  CONTINUE 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  SETUP 


SETUP  COORDINATES  OF  PANEL  NODES  AND  SLOPES  OF  PANELS 
COORDINATES  ARE  READ  FROM  INPUT  DATA  FILE  UNLESS 

THE  AIRFOIL  IS  OF  NACA  XXXX  OR  NACA  230XX  TYPE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


NACA(I) 

TAU(I) 

EPSMAX(I) 

PTMAX(I) 


DUMMY  NACA  NUMBER  FOR  TRANSFER  TO  SUBROUTINES 
BODY  AND  NACA 

DUMMY  MAX  THICKNESS  FOR  TRANSFER  TO  SUBROUTINES 
BODY  AND  NACA 

DUMMY  MAX  CAMBER  FOR  TRANSFER  TO  SUBROUTINES 
BODY  AND  NACA 

DUMMY  CHORDWISE  POSN  FOR  TRANSFER  TO  SUBROUTINES 
BODY  AND  NACA 


SUBROUTINE  SETUP 

COMMON  /BOD/  IFLAG, NLOWER, NUPPER, N0DTOT,X( 202 ) ,Y( 202) , 
C0STHE(201) ,SINTHE(201) ,SS(2) ,NP1,NP2,NP3,NP4, 
NP5,XSHIFT,YSHIFT,NAIRF0,XI(202),YI(202), 
C0STHL(201) ,SINTHL(201) 

COMMON  /PARD/  NACAD(2) ,TAUD(2) ,EPSMAD(2) ,PTMAXD(2) 
COMMON  /PAR /  NACA , TAU , EPSMAX , PTMAX 


ooooo  no 


COMMON  /NUM/  PI,PI2INV 
PI  *  3. 1415926585 

PI2INV  =  .5/PI 
C 

C  SET  COORDINATES  OF  NODES  ON  BODY  SURFACE 
C 

DO  210  I-l.NAIRFO 
NODTOT  ■  NLOWER  +  NUPPER 
K*( I  - 1 )*( NODTOT+1 ) 

IF  (IFLAG  . NE.  0)  GO  TO  10 
C 

C  NACA  SERIES  AIRFOIL  CALCULATIONS 
C 

NACA-NACAD(I) 

TAU*TAUD(I) 

EPSMAX=EPSMAD( I ) 

PTMAX*PTMAXD( I ) 

NPOINT  *  NLOWER 

SIGN  =  -1.0 

NSTART  -  0 

DO  110  NSURF  ■  1,2 

DO  100  N  *  1, NPOINT 

FRACT  *  FL0AT(N-1)/FL0AT( NPOINT) 

Z  »  .5*(1.  -  COS( PI+FRACT) ) 

J  *  NSTART  +  N 

CALL  BODY(Z,SIGN,XD,YD) 

XI (J+K)  =  XD 
YI(J+K)  *  YD 
100  CONTINUE 

NPOINT  «  NUPPER 
SIGN  «  1.0 
NSTART  *  NLOWER 
110  CONTINUE 

XI (I*( NODTOT+1))  -  XI(1+K) 

120  YI(I*( NODTOT+1))  -  YI(1+K) 

GO  TO  210 
C 

C  AIRFOIL  DATA  INPUT  THAT  ARE  NOT  NACA  SERIES  AIRFOIL. 
C 

10  READ  (1,501)  (XI( J+K), J*l, NODTOT+1) 

READ  (1,501)  (YI( J+K), J-l, NODTOT+1) 

WRITE  (6,501)  (XI( J+K), J»l, NODTOT+1) 

WRITE  (6,501)  (YI( J+K), J»l, NODTOT+1) 

501  FORMAT  (6F10.6) 

210  CONTINUE 

NP1  -  NODTOT  +  1 
NP2  «  NAIRF0+NP1+1 
C 

C  SET  SLOPES  OF  PANELS  AND  CALCULATE  AIRFOIL  PERIMETER 
C 

DO  220  I  ■  l.NAIRFO 
K  -  ( I -1)*( NODTOT+1) 

SS(I)  «  0.0 

DO  200  J  ■  1, NODTOT 

DX  »  XI(J+1+K)  -  XI(J+K) 

DY  =  YKJ+l+K)  -  YKJ+R) 


DIST  «  SQRT(DX*DX  +DY*DY) 

SS(I)  «  SS(I)  +  DIST 

SINTHL( J+K)  «  DY/DIST 
COSTHL( J+K)  «  DX/DIST 
200  CONTINUE 
220  CONTINUE 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  SUBROUTINE  BODY(Z,SIGN,X,Y)  C 
C  C 
C  RETURN  COORDINATES  OF  POINT  ON  THE  BODY  SURFACE  C 
C  C 
C  Z  =  NODE -SPACING  PARAMETER  C 
C  X,Y  ■  CARTESIAN  COORDINATES  C 
C  SIGN  «  +1.  FOR  UPPER  SURFACE  C 
C  -1.  FOR  LOWER  SURFACE  C 
C  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  BODY(Z,SIGN,XD,YD) 

COMMON  /PAR/  NACA , TAU , EPSMAX , PTMAX 
IF  (SIGN  .LT.  0.0)  Z  =  1.  -  Z 
CALL  NACA45 ( Z , THICK , CAMBER , BETA) 


XD  ■  Z  -  SIGN*THICK*SIN(BETA) 

YD  *  CAMBER  +  SIGN*THICK*COS(BETA) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 

C  SUBROUTINE  NACA45(Z, THICK, CAMBER, BETA)  C 

C  C 

C  EVALUATE  THICKNESS  AND  CAMBER  C 

C  FOR  NACA  4-  OR  5-DIGIT  AIRFOIL  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  Z  .  ABSCISSA  OF  CAMBERLINE  POINT 

C  X.Y  .  LOCAL  CARTESIAN  COORDINATES  OF  NACA-AIRFOIL 

C  SIGN  .  SURFACE  INDICATOR:  SIGN  *  1,  UPPER  SURFACE 

C  —1,  LOWER  SURFACE 

C  THICK  .  THICKNESS  AT  THE  POSITION  Z 

C  CAMBER  .  CAMBER  AT  THE  POSITION  Z 

C  BETA  .  ANGLE  BETWEEN  CAMBER  LINE  AND  X-AXIS  AT  POSITION  Z 

C  TAU  .  MAXIMUM  THICKNESS  ( INPUT) 

C  EPSMAX  .  MAXIMUM  CAMBER  (INPUT) 

C  PTMAX  .  COORDINATE  POSITION  OF  MAXIMUM  CAMBER  (INPUT) 

C  . _ 


SUBROUTINE  NACA45(Z, THICK, CAMBER, BETA) 

COMMON  /PAR/  NACA, TAU, EPSMAX, PTMAX 
THICK  »  0.  0 

IF  (Z  .LT.  l.E-10)  GO  TO  100 
THICK  »  5.*TAU*(.2969*SQRT(Z)  -  Z*(.  126  +  Z*(.3537 
+  -  Z*( .  2843  -  Z*. 1015)))) 

100  IF  (EPSMAX  .EQ.  0.0)  GO  TO  130 


ooooonooonnooo  ono  ono  non  non  onn  nnn 


IF  (NACA  .GT.  9999)  GO  TO  140 

CAMBERLINE  OF  NACA  4-DIGIT  SERIES 

IF  (2  .GT.  PTMAX)  GO  TO  110 

FORWARD  PART  OF  CAMBER  LINE 

CAMBER  *  EPSMAX/ PTMAX/ PTMAX* ( 2.  *PTMAX  -  Z)*Z 
DCAMDX  «  2. *EPSMAX/PTMAX/PTMAX*( PTMAX  -  Z) 

GO  TO  120 


AFT  PART  OF  CAMBERLINE 


110 

CAMBER 

«  EPSMAX/C 1.  -PTMAX)**2*( 1.  +  Z  - 

2.  *PTMAX)*(  1, 

DCAMDX 

*  2.  *EPSMAX/( 1.  -PTMAX)**2*( PTMAX 

-  Z) 

120 

BETA 

«  ATAN( DCAMDX) 

RETURN 

130 

CAMBER 

*  0.0 

BETA 

*  0.  0 

RETURN 

CAMBERLINE  OF  NACA  5 -DIGIT  SERIES 
140  IF  (Z  .GT.  PTMAX)  GO  TO  150 

FORWARD  PART  OF  CAMBER  LINE 
W  =  Z/PTMAX 

CAMBER  »  EPSMAX*W*( (W  -  3.  )*W  +  3.  -  PTMAX) 

DCAMDX  *  EPSMAX*3. *W*( 1.  -  W) /PTMAX 
GO  TO  120 

AFT  PART  OF  CAMBERLINE 

150  CAMBER  *  EPSMAX*(1.  -  Z) 

DCAMDX  *  -  EPSMAX 

GO  TO  120 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  GAUSS(NRHS,M,NITR) 

SOLUTION  OF  LINEAR  ALGEBRAIC  SYSTEM  BY 
GAUSS  ELIMINATION  WITHOUT  PARTIAL  PIVOTING 

'A  =  COEFFICIENT  MATRIX 

NEQNS  -  NUMBER  OF  EQUATIONS 

NRHS  *  NUMBER  OF  RIGHT  HAND  SIDES 

RIGHT-HAND  SIDES  AND  SOLUTIONS  STORED  IN 
COLUMNS  NEQNS+1  THRU  NEQNS+NRHS  OF  *A 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  GAUSS(NRHS,M,NITR) 

COMMON  / BOD/  IFLAG,NLOWER,NUPPER,N0DT0T,X(202) ,Y(202) , 


ononooonooooo 


+  COSTHE(201) ,SINTHE(201) ,SS(2) ,NP1,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRF0,XI(202) , YI(202) , 

+  COSTHL(201) ,SINTHL(201) 

COMMON  /COF/  A( 201, 211) , KEQNS 
C  WRITE  (9,900)  M.NITR 

C900  FORMAT  (IX,1 M  13,'  NITR  *' ,13/) 

C  WRITE  (9,910)(I,A(I,NP3),A(I,NP4),A(I,NP5),I=1,2*N0DT0T) 

910  FORMAT  ( 1X,I5,3E14.  6) 

IF  (M.EQ.0)  KEQNS  =  NODTOT 
NEQNS  =  KEQNS  *  NAIRFO 

NP  ■  NEQNS  +  1 

NTOT  ■  NEQNS  +  NRHS 

IF  (NITR  .GT.  0)  GO  TO  160 

C 

C  GAUSS  REDUCTION 

C 

DO  150  I  -  2, NEQNS 
IM  *  I  -  1 
C 

C  ELIMINATE  (I-l)TH  UNKNOWN  FROM 

C  ITH  THRU  ( NEQNS )TH  EQUATIONS 

C 

DO  150  J  »  I, NEQNS 
R  *  A(J,IM)/A(IM,IM) 

DO  150  K  *  I, NTOT 
150  A(J,K)  *  A(J,K)  -  R*A(IM,K) 

GO  TO  170 
C 

C  GAUSSIAN  ELIMINATION  ON  ONLY  THE  RIGHT-HAND-SIDES 

C 

160  DO  180  I  *  2, NEQNS 
IM  =1-1 
DO  180  J  =  I, NEQNS 
R  «  A(J,IM)/A(IM,IM) 

DO  180  K  =  NP.NTOT 
180  A(J,K)  =  A(J,K)  -  R*A(IM,K) 

170  CONTINUE 
C 

C  BACK  SUBSTITUTION 

C 

DO  220  K  =  NP.NTOT 

A( NEQNS, K)  =  A( NEQNS, K)/A( NEQNS, NEQNS) 

DO  210  L  =  2, NEQNS 
I  =  NEQNS  +  1  -  L 

IP  =1+1 
DO  200  J  ■  IP, NEQNS 
200  A(I,K)  ■  A(I,K)  -  A(I,J)*A(J,K) 

210  A(I,K)  =  A(I,K)/A(I,I) 

220  CONTINUE 
RETURN 
END 

*  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  SUBROUTINE  VELDIS(SINALF.COSALF) 

C 

C  COMPUTE  STEADY  FLOW  PRESSURE  DISTRIBUTION  AND  VELOCITY 


105 


o  o  o  a 


o  o  o  noon  oonooooooooo 


POTENTIAL  AT  MID-POINTS  OF  PANELS  FOR  THE  STEADY  FLOW  CASE 


C  POTENTIAL  AT  MID-POINTS  OF  PANELS  FOR  THE  STEADY  FLOW  CASE 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


VTANG  . 

PHI(I)  . 

PHILE(L). . . . 
PIN(L) 
SUMC(L)  .... 


TANGENTIAL  VELO  COMP  OF  A  FIXED  POINT  OF  THE  MOVING 
LOCAL  FRAME  OF  REFERENCE. 

DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 
I-TH  PANEL 

DIFFERENCE  OF  THE  POTENTIALS  OF  THE  LEADING  EDGE  TO 
THE  LOWER  TRAILING  EDGE  FOR  THE  RESPECTIVE  AIRFOIL 
DIFFERENCE  OF  THE  POTENTIALS  AT  A  POINT  1000  CHORD 
LENGTH  UPSTREAM  OF  THE  LE  FOR  THE  RESPECTIVE  AIRFOIL 
GAMMA  ASSOCIATED  WITH  THE  INTEGRATION  OF  THE  DISTURB¬ 
ANCE  VELOCITY  AROUND  THE  WHOLE  AIRFOIL 


SUBROUTINE  VELD IS 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHE(201) ,SINTHE(201) ,SS(2) ,NP1 ,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRFO,XI(202) ,YI(202) , 

+  C0STHL(201),SINTHL(201) 

COMMON  /COF/  A( 201, 211) ,KEQNS 
COMMON  /CPD/  CP(200),SCL,T,SCM,SGAM 
COMMON  /NUM/  PI.PI2INV 

COMMON  /SING/  Q(200) ,GAMMA(2) ,QK(200) ,GAMK(2) 

COMMON  / POT/  PHI (200) ,PHIK(200) 

COMMON  /GUST/  UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON  /EXTV/  UE( 200) ,VN( 200) 

COMMON  /GEOM/  SINALF(2) ,C0SALF(2) ,0MEGA(2) ,UX(2) ,UY(2) .PIVOT, 

+  XPRM.YPRM 

DIMENSION  CKUTTA(2) ,PHITEL(2) ,PHITEU(2) 

7770  DIMENSION  CONTR2(2) 

6666  REAL  *8  PIN(2),VELX 

DIMENSION  WGHT(5),PL0C(5),SUMC(2), 

+AANP1(50,50,5) ,AANP2(50,50,5) ,BBNP1(50,50,5) , BBNP2( 50,50,5) , 
+COSTHPC 102,6), SINTHP( 102,6) , AANP4( 2) , PHIL( 102) , UGU( 100,6), 

+VGU( 100,6) ,PHILE(2) 

LOCATION  AND  WEIGHTING  VALVES  FOR  THE  GAUSSIAN  QUADRATURE  USED 
TO  INTEGRATE  FOR  THE  VELO  POTENTIAL  AROUND  THE  AIRFOIL 


DATA  WGHT/.  11846344,.  23931434,. 28444444, 

+  .23931434,. 11846344/ 

DATA  PLOC/.  04691008,. 23076535,-50000000, 

+  .76923466,. 95300899/ 

FIND  VT  AND  CP  AT  MID-POINT  OF  I-TH  PANEL 

DO  140  L«  l.NAIRFO 

KI  -  (L-1)*NP1 

LI  -  (L-1)*N0DT0T 

SUMC(L)  -  0.0 

DO  130  I  ■  l.NODTOT 

XMID  ■  .  5*(X(I+KI)  +X( I+KI+1)) 

YMID  -  .  5*( Y( I+KI )  +Y( I+KI+1)) 

DX  «  X( I+KI+ 1 ) -X( I+KI ) 

DY  «  Y( I+KI+1) -Y( I+KI) 
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DIST  »  SQRT( DX*DX+DY*DY ) 

VTANG  =  COSALF( L)*COSTHL( I+KI )  +  SINALF(L)*SINTHL( I+KI) 

VNORM  =  SINALF( L)*COSTHL( I+KI )  -  COSALF(L)*SINTHL(I+KI) 
VTFREE  =  VTANG 

VACT  *  VTANG 

ADD  CONTRIBUTION  OF  J-TH  PANEL 

1.  CONTRIBUTION  OF  J-TH  PANEL  TO  THE  VELO  COMP  OF  THE  MIDPT 
OF  THE  I-TH  PANEL 

2.  CONTRIBUTION  TO  THE  VELO  POTENTIAL.  THIS  IS  DONE  BY 
INTEGRATING  OVER  SMALLER  PANELS  OF  THE  AIRFOIL. 


100 


120 

ISO 


155 


130 

140 


DO  155  K  »  1,5 

DX  *  PLOC(K)*(X( I+KI+1)-X( I+KI)) 

DY  *  PLOC(K)*(Y(I+KI+l)-Y(I+KI)) 

XMID  «  X( I+KI )  +  DX 

YMID  =  Y(I+KI)  +  DY 

VDUM  =0.0 
DO  150  LM  =  l.NAIRFO 
KJ  =  (LM-1)*NP1 

LJ  =  (LM-1)*N0DT0T 

DO  120  J  =  l.NODTOT 

FLOG  =0.0 
FT  AN  ■  PI 

IF  (J+KJ  .EQ.  I+KI)  GO  TO  100 
DXJ  =  XMID  -  X( J+KJ) 

DXJP  *  XMID  -  X( J+KJ+1) 

DYJ  *  YMID  -  Y( J+KJ) 

DYJP  =  YMID  -  Y( J+KJ+1) 

FLOG  «  .5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

FTAN  «  ATAN2( DY JP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE( I+KI )*C0STHE( J+KJ)  +  SINTHEC I+KI )*SINTHE( J+KJ) 
STIMTJ  »  SINTHE( I+KI )*C0STHE( J+KJ)  -  COSTHEC I+KI )*SINTHE( J+KJ) 
AA  =  PI2INV*(FTAN*CTIMTJ  +  FLOG* STIMTJ) 

B  »  PI2INV*(FLOG*CTIMTJ  -  FTAN+STIMTJ) 

VDUM  «  VDUM  -  B*Q( J+LJ)  +GAMMA(LM)*AA 

IF(K  .EQ.  3)  VACT  «  VACT  -  B*Q(J+U)  +  GAMMA(LM)*AA 

IF(K  .EQ.  3)  VNORM  *  VNORM  +  AA*Q(J+U)  +  GAMMA(LM)*B 

CONTINUE 

CONTINUE 

VTANG  «  VTANG  +VDUM  *WGHT(K) 

SUMC(L)  =  SUMC(L)  +VDUM*DIST*WGHT(K) 

CONTINUE 

PHKI+LI)  =  (  VTANG- VTFREE  )*DI  ST 

CP(I+LI)»  1.0  -  VACT+VACT 

UE(I+LI)»  VACT 

VN(I+LI)»  VNORM 

CONTINUE 

CONTINUE 


COMPUTE  DISTURBANCE  POTENTIAL  AT  THE  LEADING  EDGE  BY  LINE 
INTEGRAL  OF  THE  VELOCITY  FIELD 
FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 

DO  55  L  =  l.NAIRFO 
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YMID  =  PIVOT  *  SINALF(L)  +  (L-1)*YSHIFT 
XLE  =  -P I VOT*COSALF ( L ) +(  L- 1 ) *XSHIFT 
XL  =  XLE 

C  XL  »  0.0 

NPHI  ■  10  *  NLOWER 

PIN(L)  =0.0 

DO  30  I  =  l.NPHI 

FRACT  =  FLOAT ( I ) /FLOAT( NPHI ) 

XLP  =  -100.0  *  (1.0  -  C0S(0.5*PI*FRACT))+XLE 

IF  (I  .EQ.  1)  XLP  «  -0. 000197+XLE 
C  XLP  »  -10.0  *  (1.0  -  COS(  0.  5*PI*FRACT) ) 

DELX  =  XL  -  XLP 

XMID  -  0. 5*(XL+XLP) 

C  XMID  =  0. 5*(XL+XLP)*COSALF(L) 

C  YMID  =  0. 5*(XL+XLP)*SINALF(L) 

XL  *  XLP 

VELX  =  UGUST 

C 

C  ADD  CONTRIBUTION  OF  J-TH  PANEL 
C 

DO  40  LN  ■  l.NAIRFO 
KJ  «  (LN-1)*NP1 

U  =  ( LN- 1 )*NODTOT 

DO  20  J  »  l.NODTOT 

DXJ  ■  XMID  -  X( J+KJ) 

DXJP  =  XMID  -  X( J+KJ+1) 

DYJ  *  YMID  -  Y( J+KJ) 

DYJP  »  YMID  -  Y( J+KJ+1) 

FLOG  *  .  5*ALOG( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 

FT AN  -  ATAN2( DY JP*DXJ-DXJP*DYJ , DXJP*DXJ+DY JP*DY J) 

CALMTJ  -  -COSALF(L)*COSTHE( J+KJ)  -  SINALF(L)*SINTHE( J+KJ) 

SALMTJ  =  -SINALF( L)*COSTHE( J+KJ)  +  COSALF( L)*SINTHE( J+KJ) 

APY  ■  PI2INV*(FTAN*CALMTJ  +  FLOG* SALMTJ) 

BPY  -  PI2INV*(FLOG*CALMTJ  -  FTAN*SALMTJ) 

VELX  =  VELX  -  DPROD(BPY,Q( J+LJ))  +DPROD(GAMMA(LN) ,APY) 

20  CONTINUE 

40  CONTINUE 

PIN(L)  «  PIN(L)  +  VELX  *  DBLE(DELX) 

30  CONTINUE 

55  CONTINUE 

C 

C  COMPUTATION  OF  THE  VELOCITY  POTENTIAL  FOR  MIDPOINT  OF  EACH  PANEL 
C 

DO  240  L  *  l.NAIRFO 

LI  =  (L-1)*NODTOT 

PHP  «  -PIN(L) 

C 

C  BEGIN  WITH  LOVER  SURFACE 
C 

DO  230  I  »  NLOWER, 1,-1 

PHC  ■  PHP-PHKI+LI) 

PHKI+LI)  -  0. 5*(PHP+PHC) 

230  PHP  -  PHC 

PHITEL(L)  «  PHC 

C 


C  RESET  FOR  UPPER  SURFACE 
C 


PHP  =  -PIN(L) 

DO  250  I  =  NLOWER+1 .NODTOT 
PHC  ■  PHP+PHKI+LI) 

PHI(I+LI)  *  0. 5*(PHP+PHC) 

250  PHP  «  PHC 

PHITEU(L)  *  PHC 
240  CONTINUE 

DO  334  L  ■  1  ,  NAIRFO 
WRITE  (6,1010)  L 
WRITE  (6,1000) 

KI  «  (L-l)  *  NP1 

LI  »  (L-l)  *  NODTOT 

SUMC(L)  «  SUMC(L)/SS(L) 

DO  333  I  «  1  ,  NODTOT 

XMID  «  .  5*(X( I+KI)  +  X( I+KI+1)) 

YMID  *  . 5*(Y(I+KI)  +  Y( I+KI+1)) 

333  WRITE  (6,1050)  I+LI , XMID, YMID, Q( I+LI ) ,GAMMA(L) ,CP( I+LI) ,UE( I+LI) , 

+  VN( I+LI), PHI (1+LI),SUMC(L) 

334  WRITE  (6,235)  PIN(L) 

235  FORMAT  (lX.’PHI  AT  LEADING  EDGE  ■' ,F10.  6,/) 

1000  FORMAT( / ,4X, 1 J' ,4X,'X(J)' ,5X,'Y(J)f ,6X,’Q(J)' ,6X, ’GAMMA’ ,5X, 

+  'CP( J)' ,6X, ’ V( J) ' ,6X, ' VN( J)' ,6X, 'PHI ' ,4X, ' INTGAMMA' ,/) 

1010  FORMAT(//'  AIRFOIL  NUMBER' ,14,/) 

1050  FORMAT( 15 ,9F10.  6) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c  SUBROUTINE  FANDM( SINALF.COSALF) 

C 

C  COMPUTE  AND  PRINT  OUT  CD, CL, CM 

C  INTEGRATE  PRESSURE  DISTRIBUTION  BY  TRAPEZOIDAL  RULE 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c  CP(I)  .  PRESSURE  COEFFICIENT  OF  THE  I-TH  PANEL 

C  CL(L)  .  COEFFICIENT  OF  LIFT  FOR  THE  L-TH  AIRFOIL 

C  CD(L)  .  COEFFICIENT  OF  DRAG  FOR  THE  L-TH  AIRFOIL 

C  CM(L)  .  COEFFICIENT  OF  MOMENT  FOR  THE  L-TH  AIRFOIL  WITH 

C  RESPECT  TO  THE  LEADING  EDGE 

C  CFX(L)  .  COEFFICEIENT  OF  TOTAL  FORCE  IN  X-DIR  OF  GLOBAL  SYS. 

C  CF¥(L)  .  COEFFICEIENT  OF  TOTAL  FORCE  IN  Y-DIR  OF  GLOBAL  SYS. 

C  _ 

c 


SUBROUTINE  FANDM(SIN1,SIN2,C0S1,C0S2) 

COMMON  /BOD/  IFLAG, NLOWER.NUPPER, NODTOT, X(202) ,Y( 202) , 

+  C0STHE(201) ,SINTHE(201) ,SS(2) ,NP1,NP2,NP3,NP4, 

+  NP5, XSHIFT, YSHIFT, NAIRFO, XI ( 202 ),YI( 202), 

+  COSTHL(201) ,SINTHL(201) 

COMMON  /CORV/  CV(2,200) ,XC(2,200) ,YC(2,200) ,M,TD,CWX(2,200) , 
+  CWY(2,200),XCI(  2,200)  ,YCI(2, 200) 

COMMON  /CPD/  CP(200) ,SCL,T,SCM,SGAM 
COMMON  /8ING/  Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 

COMMON  /GEOM/  SINALF(2) ,C0SALF(2) ,0MEGA(2) ,UX(2) ,UY(2) , PIVOT, 
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+  XPRM , YPRM 

DIMENSION  CM( 2) ,CD(2) ,CL(2) ,CFX(2) ,CFY(2) 

INITIALISE  COEFFICIENT 

DO  110  L  =  l.NAIRFO 
CM(L)  =0.0 
CFX(L)  =0.0 
110  CFY(L)  =0.0 

INTEGRATE  AROUND  THE  AIRFOIL  TO  GET  GLOBAL  X  AND  Y  FORCES 

DO  120  L  =  l.NAIRFO 
KI  =  (L-1)*(N0DT0T+1) 

LI  =  (L-l)*NODTOT 

DO  100  I  =  l.NODTOT 

XMID  «  .  5*(XI( I+KI)  +  XI(I+KI+1)) 

YMID  =  .  5*(YI( I+KI )  +  YKI+KI+l)) 

DX  =  XI( I+KI+1)  -  XI ( I+KI) 

DY  =  YIC I+KI+1)  -  YKI+KI) 

CFX(L)  =  CFX(L)  +  CP(I+LI)*DY 
CFY(L)  *  CFY(L)  -  CP(I+LI)*DX 
CM(L)  =  CM(L)  +  CP( I+LI )*( DX+XMID  +  DY+YMID) 

100  CONTINUE 
120  CONTINUE 

DECOMPOSE  INTO  LIFT  AND  DRAG  COMPONENTS  W.R.T.  RESP  LOCAL  SYSTEM 
DO  130  L  =1 .NAIRFO 

CD(L)  =  CFX( L)*COSALF( L)  +  CFY(L)*SINALF(L) 

CL(L)  »  CFY(L)*COSALF(L)  -  CFX(L)*SINALF(L) 

130  WRITE  (6,1000)  L,CD(L) ,CL(L) ,CM(L) 

IF  (M  .EQ.  0)  RETURN 
IF  (SCL  .NE.  0.0)  THEN 

WRITE  (8,1100)  T,CL( 1)/SCL,CL( 2)/SCL,CM( 1)/SCM,CM(2)/SCM,CD( 1) , 
+CD( 2) ,GAMK( 1)/SGAM,GAMK( 2)/SGAM 
ELSE 

WRITE  (8,1100)  T,CL( 1) ,CL( 2) ,CM( 1) ,CM(2) ,CD( 1) ,CD( 2) 

ENDIF 

1000  FORMAT(// , '  AIRFOIL  NO  ',14./,'  CD  »’,F10.  6, 

+  1  CL=' ,F10. 6 , '  CM*' ,F10.6) 

1100  FORMAT(9F10.  6) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  I NFL  (NITR) 

CALCULATE  INFLUENCE  COEFFICIENTS 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

INFLUENCE  COEFFICIENTS  ON  THE  AIRFOIL  DUE  TO  THE  AIRFOIL: 

AAN(I.J)  _  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I -TH  PANEL  DUE 

TO  A  SOURCE-DIST  OF  UNIT  STRENGTH  ON  THE  J-TH  PANEL 
SUMAAN(I.L)  . .  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
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SOURCE  DIST  OF  UNIT  STRENGTH  FROM  THE  L  AIRFOIL 

BBN(I,J)  -  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 

TO  A  VORTEX-DIST  OF  UNIT  STRENGTH  ON  THE  J-TH  PANEL 

SUMAAN(I,L)  ..  NORMAL  VELO  AT  rHE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
VORTEX  DIST  OF  UNIT  STRENGTH  FROM  THE  L  AIRFOIL 

INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT: 

AYNPl(L.J)  ...  Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRE¬ 
NGTH  FROM  THE  J-TH  PANEL 

AYNP1(L,NP3)  .  Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  PANEL  DUE  TO  A  SOURCE  DIST  OF  UNIT  STREN¬ 
GTH  FROM  THE  WAKE  PANEL  OF  THE  OTHER  AIRFOIL. 

(USED  ONLY  FOR  BXNP1(L,NP3)  SINCE  THERE  IS  NO  SOURC 
DIST  ON  THE  WAKE  PANEL) 

BYNPl(L.J)  ...  Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  VORTEX  DIST  OF  UNIT  STRE¬ 
NGTH  FROM  THE  J-TH  PANEL 

BYNP1(L,NP3)  .  Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  A  VORTEX  DIST  OF  UNIT  STR¬ 
ENGTH  FROM  THE  WAKE  PANEL  OF  THE  OTHER  AIRFOIL. 

CYNPl(L.N)  . . .  Y-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  THE  N-TH  CORE  VORTEX  OF 
UNIT  STRENGTH 

CXNPl(L.N)  . . .  X-VELO  COMP  AT  THE  MIDPOINT  OF  THE  WAKE  PANEL  FROM 
THE  L-TH  AIRFOIL  DUE  TO  THE  N-TH  CORE  VORTEX  OF 
UNIT  STRENGTH 

INFLUENCE  COEFFICIENTS  ON  THE  AIRFOIL  DUE  TO  THE  WAKE 

SUMCCN(I)  ...  NORMAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL  DUE 
TO  ALL  POINT  VORTICES  OF  ACTUAL  STRENGTH. 

SUMCCT(I)  ...  TANGENTIAL  VELO  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL 
DUE  TO  ALL  POINT  VORTICES  OF  ACTUAL  STRENGTH. 


SUBROUTINE  INFL  (NITR) 

COMMON  / BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRF0,XI(202) ,YI(202) , 

+  COSTHL(201),SINTHL(201) 

COMMON  /NUM/  P1.PI2INV 

COMMON  /WAK/  VYW(2) ,VXW(2) ,WAKE(2) ,DT 

COMMON  /CORV/  CV(2,200) ,XC(2,200) ,YC(2,200) ,M,TD,CWX(2,200) , 
CWY(2 ,200)  ,XCI(2,200)  ,YCI(2,200) 

COMMON  /INFl/  AAN(201,201) ,BBN(201,201) ,AYNP1(2,201) ,BYNP1(2,201) , 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201) ,SUMCCT(201) ,CYNP1(2,200) , 

+  CXNP1(2,200) 

COMMON  /GEOM/  SINALF(2) ,COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2) .PIVOT, 

+  XPRM.YPRM 

DIMENSION  GAYNP1(2,201) ,GBYNP1(2,201) 

IF  (M  .GT.  1)  GO  TO  510 

INFLUENCE  COEFFICIENT  ON  THE  I-TH  PANEL  BY  THE  J-TH  PANEL  FROM 
THE  SAME  AIRFOIL. 


Ill 


noon 


DO  200  L  =  1 ,NAIRF0 

LI  =  (L-1)*N0DT0T 

KI  =  (L-1)*NP1 

DO  120  I  =  1.N0DT0T 

XMID  -  .  5*(XI( I+KI)  +  XICI+KI+D) 

YMID  =  . 5*( YI( I+KI)  +  YI( I+KI+1)) 

SUMAAN(I+LI,L)  =0.0 
SUMBBN(I+LI,L)  =0.0 
DO  110  J  =  1.N0DT0T 
FLOG  =  0. 0 

FT AN  =  PI 

IF  (J+LI  .EQ.  I+LI)  GO  TO  100 
DXJ  =  XMID  -  XI (J+KI) 

DXJP  =  XMID  -  XKJ+KI+l) 

DYJ  =  YMID  -  YI(J+KI) 

DYJP  =  YMID  -  YI( J+KI+1) 

FLOG  =  .  5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 

FT AN  =  AT AN  2 ( DYJP*DXJ - DXJP*DYJ , DX JP*DX J+DY JP*DY J ) 

100  CTIMTJ  =  COSTHL( I+KI )*C0STHL( J+KI )  +  SINTHL(I+KI)*SINTHL( J+KI) 
STIMTJ  =  SINTHL( I+KI )*COSTHL( J+KI)  -  C0STHL( I+KI )*SINTHL( J+KI) 
AAN( I+LI, J+LI)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG+STIMTJ) 

BBN( I+LI, J+LI)  =  PI2INV*(FLOG*CTIMTJ  -  FTAN* STIMTJ) 

SUMAAN( I+LI ,L)  =  SUMAAN( I+LI ,L)  +  AAN( I+LI , J+LI) 

SUMBBN( I+LI , L)  =  SUMBBN( I+LI , L)  +  BBN( I+LI , J+LI) 

110  CONTINUE 
120  CONTINUE 
200  CONTINUE 
510  CONTINUE 

IF  (NITR  .GT.  0)  GO  TO  271 

INFLUENCE  COEFFICIENT  ON  THE  I-TH  PANEL  BY  THE  J-TH  PANEL  FROM 
THE  OTHER  AIRFOIL. 

DO  270  L  =  l.NAIRFO 
LI  =  (L-l)*NODTOT 

KI  =  (L-1)*NP1 

DO  260  I  =  l.NODTOT 
XMID  ■  .  5*(X( I+KI)  +  X( I+KI+1)) 

YMID  =  . 5*(Y(I+KI)  +  Y( I+KI+1)) 

SUMAAN( I+LI , 3-L)  =0.0 
SUMBBN( I+LI , 3-L)  =0.0 
DO  250  J  «  NODTOT+2 , 2+NODTOT+l 
DXJ  =  XMID  -  X(J-KI) 

DXJP  =  XMID  -  X( J-KI+1) 

DYJ  «  YMID  -  Y(J-KI) 

DYJP  =  YMID  -  Y( J-KI+1) 

FLOG  ■  .  5*AL0G( ( DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 

FTAN  »  ATAN2 ( DY JP+DXJ -DXJP+DY J , DXJP+DXJ+DY JP+DY J ) 

CTIMTJ  ■  COSTHE( I+KI )*COSTHE( J-KI )  +  SINIHE(I+KI)*SINIHE( J-KI) 
STIMTJ  ■  SINTHE( I+KI )*COSTHE( J-KI)  -  COSTHE( I+KI )*SINTHE( J-KI) 
AAN( I+LI, J- LI -1)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG* STIMTJ) 

BBN( I+LI , J-LI - 1 )  =  PI2INV*( FLOG* CTIMTJ  -  FTAN* STIMTJ) 

SUMAANC I+LI, 3-L)  =  SUMAAN( I+LI, 3-L)  +  AAN(I+LI,J-LI-1) 

SUMBBN( I+LI, 3-L)  ■  SUMBBN( I+LI, 3-L)  +  BBN(I+LI,J-LI-1) 

250  CONTINUE 
260  CONTINUE 
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270  CONTINUE 

END  OF  STEADY  FLOW  CALAULATION  FOR  INFLUENCE  COEFFICIENT. 

IF  (M.EQ.O)  RETURN 

271  CONTINUE 

INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT  FROM  THE  AIRFOIL. 

DO  130  L  «  l.NAIRFO 

I  «  NP1*L 

XMID  *  .5*(X(I)  +  X(NAIRF0*NP1+L)) 

YMID  =  . 5*(Y( I)  +  Y(NAIRF0*NP1+L)) 

DO  130  MM  »  l.NAIRFO 
LJ  *=(MM-1)*N0DT0T 

KJ  «(MM-1)*NP1 

DO  130  J  *  l.NODTOT 
DXJ  «  XMID  -  XC J+KJ) 

DXJP  *  XMID  -  XC J+KJ+1) 

DYJ  *  YMID  -  Y( J+KJ) 

DYJP  *  YMID  -  Y( J+KJ+1) 

FLOG  =  .  5*AL0G( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 

FT AN  =  ATAN2(DYJP*DXJ“DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  *  C0STHE( I )*C0STHE( J+KJ )  +  SINTHEC I )*SINTHE( J+KJ) 

STIMTJ  «  SINTHE( I )*COSTHE( J+KJ)  -  COSTHEC I)*SINTHE( J+KJ) 

GAYNP 1 ( L , J+LJ )  =  P I 2 I NV* ( FTAN*COSTHE ( J+KJ )  -  FLOG*SINTHE( J+KJ) ) 
GBYNP1(L,J+LJ)  =  P I 2 I NV* ( FLOG+COSTHE ( J+KJ )  +  FTAN*SINTHE( J+KJ) ) 
AYNP1(L,J+LJ)  »  GAYNP1  ( L , J+LJ)*COSALF( L) -GBYNP1  ( L, J+U)*SINALF( L) 
BYNP1(L,J+U)  =  GAYNP1(L,J+U)*SINALF(L)+GBYNP1(L,J+U)*C0SALF(L) 
130  CONTINUE 

INFLUENCE  COEFICIENT  OF  WAKE  ELEMENT  DUE  TO  WAKE  ELEMENT  FROM 
THE  OTHER  AIRFOIL 

DO  300  L  *  l.NAIRFO 

I  =  NP1*L 

XMID  =  . 5*(X( I)  +  X(NAIRF0*NP1+L)) 

YMID  «  . 5*(Y( I)  +  Y(NAIRF0*NP1+L)) 

KJ  =  (L-1)*NP1 

MJ  =  2*NP1 

NJ  *  MJ+3 

DXJ  «  XMID  -  X(MJ-KJ) 

DXJP  *  XMID  -  X(NJ-L) 

DYJ  »  YMID  -  Y(MJ-KJ) 

DYJP  *  YMID  -  Y(NJ-L) 

FLOG  *  .  5*AL0G( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ) ) 

FTAN  *  ATAN2 ( DYJP+DXJ - DXJP+DYJ , DX  JP*DXJ+DY JP*DYJ ) 

CTIMTJ  -  COSTHE ( I ) *COSTHE ( M J -KJ )  +  SINTHE( I)*SINTHE(MJ-KJ) 

STIMTJ  *  SINTHEC I )*COSTHE(MJ-KJ)  -  COSTHEC I) *SINTHE(MJ-KJ) 
GAYNP1CL.NP3)  «  PI2INV*CFTAN*C0STHE(MJ-KJ)  -  FLOG*SINTHE(MJ-KJ)) 
GBYNP1CL.NP3)  ■  PI2INV*(FLOG*COSTHE(MJ-KJ)  +  FTAN*SINTHECMJ-KJ) ) 
AYNP1CL.NP3)  *  GAYNP1CL,NP3)*C0SALFCL)-GBYNP1CL,NP3)*SINALFCL) 
BYNP1CL.NP3)  ■  GAYNP1CL,NP3)*SINALFCL)+GBYNP1CL,NP3)*C0SALFCL) 
300  CONTINUE 
C 
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C  INFLUENCE  COEFFICIENT  ON  THE  AIRFOIL  BY  THE  WAKE  ELEMENT. 

C 

DO  160  L  =  l.NAIRFO 

LI*( L- 1 )*NODTOT 

KI*(L-1)*NP1 

DO  140  I  »  l.NODTOT 

XMID  «  . 5*(X(I+KI)  +  X( I+KI+1)) 

YMID  »  .  5*(Y( I+KI)  +  Y( I+KI+1)) 

DO  145  MM  a  i,  NAIRFO 
J  «  NP1*MM 

DXJ  «  XMID  -  X(J) 

DXJP  *  XMID  -  X(NAIRF0*NP1+MM) 

DYJ  *  YMID  -  Y(J) 

DYJP  *  YMID  -  Y(NAIRF0*NP1+MM) 

FLOG  a  .  5*AL0G( ( DX JP*DX JP+DY JP*DY JP ) / ( DXJ*DX J+DY J*DY J ) ) 

FT AN  a  AT AN 2 C DY JP*DXJ -DXJP*DY J , DXJP*DXJ+DY JP*DY J ) 

CTIMTJ  «  COSTHE( I+KI)*COSTHE( J)  +  SINTHE(I+KI)*SINTHE( J) 
STIMTJ  a  SINTHE( I+KI)*COSTHE( J)  -  COSTHE( I+KI)*SINTHE( J) 

AAN( I+LI , 2*N0DT0T+MM)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
BBN(I+LI,2*N0DT0T+MM)  =  PI2INV*(FL0G*CTIMTJ  -  FTAN*STIMTJ) 

145  CONTINUE 
140  CONTINUE 
160  CONTINUE 
C 

IF  (M.EQ.  1)  RETURN 
MM1  =  M  -  1 
C 

C  INFLUENCE  COEFFICIENT  ON  THE  WAKE  ELEMENT  BY  THE  CORE  VORTICES. 
C 

DO  350  L  »  l.NAIRFO 

XMID  »  0.  5*(X(NP1*L)  +  X(NP1*NAIRF0+L)) 

YMID  a  0.  5*(Y(NP1*L)  +  Y(NP1*NAIRF0+L) ) 

DO  240  MM  =  1,  NAIRFO 
KN  »  (MM-1)*MM1 
DO  230  N  =  1,MM1 
DX  =  XMID  -  XC(MM,N) 

DY  »  YMID  -  YC(MM,N) 

DIST2 ,  =  DX*DX+DY*DY 

GCYNP1  «  -PI2INV*DX/DIST2 
GCXNP1  a  +P I 2 INV*DY / D I ST2 

CYNP1(L,N+KN)  *  GCYNP1*C0SALF(L)+GCXNP1*SINALF(L) 

CXNP1(L,N+KN)  =  -GCYNP1*SINALF(L)+GCXNP1*C0SALF(L) 

230  CONTINUE 
240  CONTINUE 
350  CONTINUE 

IF  (NITR.GT.  0)  RETURN 
C 

C  INFLUENCE  COEFFICIENT  ON  THE  AIRFOIL  BY  THE  CORE  VORTICES 
C 

DO  400  L>  l.NAIRFO 

LI  «(L-1)*N0DT0T 

KI  «(L-1)*NP1 

DO  220  I  »  l.NODTOT 

XMID  »  0.  5*(X(I+KI)  ♦  X( I+KI+1)) 

YMID  =  0. 5*(Y(I+KI)  +  Y( I+KI+1)) 


oooooooooononoooonnoooooo 


SUMCCN(I+LI)  «  0.0 
SUMCCT( I+LI )  -0.0 
DO  280  MM  »  1,  NAIRF0 
KN  -  (MM-1)*MM1 
DO  210  N  *  1,MM1 
DX  -  XMID  -  XC(MM,N) 

DY  *  YMID  -  YC(MM,N) 

DIST  -  SQRT(DX*DX+DY*DY) 

COSTHN  «  DX/DIST 

SINTHN  ■  DY/DIST 

CTIMTN  ■  COSTHEC I+KI )*C0STHN  +  SINTHE(I+KI)*SINTHN 

STIMTN  ■  SINTHE( I+KI )*C0STHN  -  COSTHEC  H*I)*SINTHN 

CCN  «  -CTIMTN/ DIST 
CCT  *  -STIMTN/DIST 

SUMCCNC I+LI )  ■  SUMCCNC I+LI)  +  CCN*CV(MM,N) 

SUMCCT( I+LI )  «  SUMCCT( I+LI)  +  CCT*CV(MM,N) 

210  CONTINUE 
280  CONTINUE 

SUMCCN( I+LI)  *  PI2INV*SUMCCN( I+LI ) 

SUMCCT( I+LI)  =  PI2INV*SUMCCT( I+LI ) 

220  CONTINUE 
400  CONTINUE 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  COEF  ( SINALF, COS ALF, OMEGA, UX,UY,NITR) 

SET  COEFFICIENTS  OF  N  EQUS  ARISING  FROM  FLOW 

TANGENCY  CONDITIONS  AT  MID  POINTS  OF  PANELS 
SOLVING  THE  N- SOURCE  STRENGTHS  IN  TERMS  OF  THE 
VORTICITY  STRENGTH  (RESULTING  IN  2  RHS) 

KUTTA  CONDITION  IS  SATISFIED  SEPARATELY  TO  OBTAIN 
THE  VORTICITY  STRENGTH 

THIS  SOLUTION  METHOD  IS  DESIRED  FOR  UNSTEADY  FLOW 


A( I , J)  .  LHS,  NORMAL  VELOCITY  AT  THE  MIDPOINT  OF  THE  I-TH 

PANEL  INDUCED  BY  UNIT  SOURCE  DIST  ON  THE  J-TH  PANEL 

A( I ,NP3)  -  FIRST  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

PANEL  WHICH  IS  DEPENDENT  ON  VORTICITY  STRENGTH  OF 
THE  FIRST  AIRFOIL 

A(I,NP4)  -  SECOND  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

PANEL  WHICH  IS  DEPENDENT  ON  VORTICITY  STRENGTH  OF 
THE  SECOND  AIRFOIL 

A( I ,NP5)  -  THIRD  RHS,  COMPONENT  OF  NORMAL  VELO  AT  THE  I-TH 

PANEL  WHICH  IS  INDEPENDENT  OF  THE  CIRCULATION  OF 
BOTH  AIRFOILS 


SUBROUTINE  COEF  (NITR) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHEC 201), SINTHE( 201 ),SS( 2), NP1,NP2,NP3,NP4, 

+  NP5,XSHIFT,YSHIFT,NAIRF0,XI(202) ,YI(202) , 

COSTHL(201),SINTHL(201) 
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mammmi 


oonnnnoonon 


non  onn  non  n  o  o 


COMMON  /COF/  A( 201, 211) ,KEQNS 

COMMON  /SING/  Q(200) ,GAMMA(2) ,QK(200) ,GAMK(2) 

COMMON  /WAX/  VYW(2)  ,VXW(2)  ,WAKE(2)  ,DT 

COMMON  /CORV/  CV(2,200)  ,XC(2,200)  ,YC(2,200)  ,M,TD,CWX(2,200)  , 

+  CWY(2,200)  ,XCI(2,200)  ,YCI(2,200) 

COMMON  /INF1/  AAN(201 ,201) JBBN(201,201) ,AYNP1(2,201) ,BYNP1(2,201) , 
+  SUMAAN(201,2) , SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201),SUMCCT(201),CYNP1(2,200), 

+  CXNP1(2 ,200) 

COMMON  /GUST/  UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON  /MATRIX/  AMAT(201,211) 

COMMON  /GEOM/  SINALF(2) ,C0SALF(2) ,OMEGA(2) ,UX(2) ,UY(2) , PIVOT, 

+  XPRM.YPRM 

NEQS  «  NODTOT*  NAIRFO 
IF  (NITR  .  GT.  0)  GO  TO  91 

SET/RESET  LHS  MATRIX  A(I,J)  FOR  EACH  NEW  TIME  STEP 

DO  110  I  «  1 , 2+NODTOT 
DO  110  J  ■  1 , 2*N0DT0T 
110  A(I,J)  *  AAN( I , J) 

IF  (M.NE. 0)  GOTO  91 

FILL  IN  THE  RIGHT  HAND  SIDE  FOR  STEADY  FLOW 

DO  310  L=l, NAIRFO 

LI  =  (L-l)*NODTOT 

KI  =  (L-1)*NP1 

DO  310  I  »  1,  NODTOT 

XMID  -  0.5  *  (XKI+KI)  +  XIU+KI+1)) 

YMID  ■  0.5  *  (YICI+KI)  +  YKI+KI+l)) 

A( I+LI ,NP3)  *  -SUMBBN(I+LI,1) 

A( I+LI ,NP4)  «  -SUMBBN( I+LI ,2) 

310  A( I+LI ,NP5)  *  SINTHLC I+KI)*COSALF(L)  -  COSTHL(I+KI)*SINALF(L) 
RETURN 


FILL  IN  THE  RIGHT  HAND  SIDE  FOR  UNSTEADY  FLOW 

91  DO  219  L*l, NAIRFO 

LI  =  (L-l)*NODTOT 

KI  «  (L-1)*NP1 

DO  210  I  «  1, NODTOT 

XMID  -  0.5  *  (XKI+KI)  +  XKI+KI+l)) 

YMID  *  0.5  *  (YKI+KI)  +  YId+KI+1)) 

A( I+LI.NP3)  »  -SUMBBNC I+LI , 1 )  +  BBN(I+LI,NP3)*SS(1)/VAKE(1) 

A( I+LI ,NP4)  «  -SUMBBNC I+LI , 2 )  +  BBN(I+LI,NP4)*SS(2)/VAKE(2) 

A( I+LI.NP5)  ■  -BBN( I+LI , NP3 )+GAMMA( 1 )*SS( 1 ) /WAKE( 1 ) -BBN( I+LI ,NP4 ) 
+*GAMMA(2)*SS(2)/WAKE(2)  +  SINTHLC I+KI )*  ((l.-KJG(I+LI))*COSALF(L)- 
+VG( I+LI )*SINALF( L)+UX( L) )  -  COSTHL(I+KI)*  ((1. +UG(I+LI)) 

+*SINALF( L)+VG( I+LI)*COSALF( L)+UY( L) )+OMEGA( L)*(YMID*SINTHL( I+KI ) 

+  +  XMID*COSTHL( I+KI ) ) 

ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .EQ.  1)  GOTO  210 

A( I+LI ,NP5)  «  ACI+LI.NP5)  -  SUMCCN(I+LI) 


no  ooononoonooooonnooonoooo 


210  CONTINUE 

DO  300  1*1 , 2*N0DT0T 
DO  300  J*1 , 2*N0DT0T+3 
300  AMAT(I,J)  *  A( I , J) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  KUTTA  ( ALPHA, SINALF.COSALF, OMEGA, UX,UY) 

USING  KUTTA  CONDITION  TO  DETERMINE  VORTICITY 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


GAMMA 


B1(I)  •• 
B2( I)  .. 
B3(I) 
AA(I,J) 


BBC  I)  .. 


CIRCULATION  FOR  THE  STEADY  FLOW  CASE  AND  ALSO 
CIRCULATION  FOR  THE  PREVIOUS  TIME  STEP  FOR  UNSTEADY  CASE 
CIRCULATION  AT  THE  CURRENT  TIME  STEP  FOR  UNSTEADY  CASE 
SOURCE  STRENGTH  AT  THE  CURRENT  TIME  STEP 
QK(I)  «  B1(I)*GAMK(1)+B2(I)*GAMK(2)+B3(I) 

THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDUCED  BY  A 

CIRCULATION  OF  UNIT  STRENGTH  FROM  AIRFOIL  1 

THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDUCED  BY  A 

CIRCULATION  OF  UNIT  STRENGTH  FROM  AIRFOIL  2 

THAT  PART  OF  THE  SOURCE  STRNGTH  WHICH  IS  INDEPENDENT 

OF  THE  CIRCULATION  FROM  BOTH  AIRFOILS 

THAT  PART  OF  THE  TANGENTIAL  VELOCITY  AT  THE  TRAILING 

EDGE  PANELS  WHICH  IS  INDUCED  BY  A  CIRCULATION  OF  UNIT 

STRENGTH  BY  AIRFOIL  J 

THAT  PART  OF  THE  TANGENTIAL  VELOCITY  AT  THE  TRAILING 
EDGE  PANELS  WHICH  IS  INDEPENDENT  OF  THE  CIRCULATION 


SUBROUTINE  KUTTA(NITR.PVTAG) 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHE(201),SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRF0,XI(202) ,YI(202) , 

+  C0STHL(201) ,SINTHL(201) 

COMMON  /COF/  A( 201, 211) ,KEQNS 

COMMON  /SING/  Q(200) ,GAMMA(2),QK(200) ,GAMK(2) 

COMMON  /WAK/  VYW(2) ,VXW( 2) ,WAKE(2) ,DT 

COMMON  /CORV/  CV( 2 ,200) ,XC(2,200) ,YC(2,200) ,M,TD,CWX(2,200) , 

+  CWY(2 ,200)  ,XCI(2 ,200)  ,YCI(2,200) 

COMMON  /INF1/  AAN(201 ,201) ,BBN(201,201),AYNP1(2,201),BYNP1(2,201) , 
+  SUMAAN(201,2),SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201) ,SUMCCT(201) ,CYNP1(2,200) , 

+  CXNP1(2,200) 

COMMON  /GUST/  UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON  /MATRIX/  AMAT(201,211) 

COMMON  /BMAT/  B1,B2,B3 

COMMON  /GEOM/  SINALF(2) ,C0SALF(2) ,OMEGA(2) ,UX(2) ,UY(2) , PIVOT, 

+  XPRM.YPRM 

COMMON  /GIES/  NGIES 

DIMENSION  Bl(200) ,B2(200) ,B3(200) ,AA(4,2) ,BB(5) ,EE(2) ,FF(2) , 
+GG(2),RADI(2),AAA(2),BBB(2),CCC(2),DDD(2),E££(2),FFF(2), 

+DG( 2) ,E( 2) ,GAMA( 2) ,DELG( 2) ,DGAMK( 2 , 3) ,COEFL( 2 , 3) 


RETRIEVE  SOLUTION  FROM  A-MATRIX  FOR  THE  SOURCE  STRENGTH  IN  TERMS 
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u  u  o  u  o 


non  o  n  a  non  noon  oo 


THE  TWO  GAME  CONTRIBUTION  AND  THE  CONSTANT  COMPONENT 


COUNT  «  0 

DO  50  I  *  1,N0DT0T*2 
B1(I)  »  A(I,NP3) 

B2(I)  -  A( I ,NP4) 

50  B3(I)  *  A(I,NP5) 

IF  (M.  GT.  0)  GOTO  400 

STEADY  KUTTA  CONDITION 

COMPUTE  TANGENTIAL  VELOCITIES 
DO  425  L  «  1,NAIRF0 
LI  =  (L-1)*N0DT0T 
KI  *  (L-1)*NP1 

KK  =  (L-l)*2 
DO  430  K  *  1,2 
IF  (K  .EQ.  1)  I  -  1 

IF  (K  .EQ.  2)  I  *  NODTOT 

XMID  =  0.  5  *  (XKI+KI)  +  XKI+KI+l)) 

YMID  *  0.5  *  (YKI+KI)  +  YICI+KI+D) 

AA(K+KK,1)  *  SUMAANC I+LI , 1) 

AA(K+KK,2)  *  SUMAAN( I+LI, 2) 

BBCK+KK)  *  C0SALF(L)*C0STHL( I+KI)  +  S INALF( L)*SINTHL( I+KI ) 
DO  419  J  *  1 , 2*N0DT0T 

AA(K+KK,1)  »  AA(K+KK, 1)  -  BBN( I+LI , J)*B1( J) 

AA(K+KK,2)  *  AA( K+KK , 2 )  -  BBN(I+LI , J)*B2( J) 

BB(K+KK)  «  BBCK+KK)  -  BBNC I+LI , J)*B3C J) 

419  CONTINUE 
430  CONTINUE 
425  CONTINUE 

SET  UP  KUTTA  CONDITIONS 

DO  450  I»l,2 
LI  «  CI-D*NAIRFO+l 
DGAMKCI.l)  *  AACL1.1)  +  AACLI+1,1) 

DGAMKCI.2)  =  AACLI.2)  +  AACLI+1,2) 

450  DGAMKCI.3)  =  -CBBCLI)  +  BBCLI+1)) 

SOLVE  KUTTA  CONDITION  BY  GAUSS 

R=DGAMKC  2 , l)/DGAMKt 1 , 1) 

DO  460  K  =  2,3 

460  DGAMKC 2, K)  *  DGAMKC 2 ,K) -R+DGAMKC 1 ,K) 

GAMMAC2)  *  DGAMKC 2, 3) /DGAMKC 2, 2) 

GAMMAC1)  *  C DGAMKC 1,3) -DGAMKC 1,2)*GAMMAC 2)) /DGAMKC 1,1) 

CALCULATE  SOURCE  STRENGTH 

I  DO  470  L  =  l.NAIRFO 

T  LI  «CL-1)*N0DT0T 

DO  470  I  =  1, NODTOT 

470  QC I+LI)  -  GAMMAC1)*B1U+LI)  +  GAMMAC2)*B2(I+LI)  +  B3CI+LI) 
RETURN 


C 


non  non  n  o  o  o  n 


UNSTEADY  KUTTA  CONDITION 


FIND  VT  AT  TRAILING  EDGE  PANELS 

INITIAL  GUESS  FOR  GAMK  FROM  PREVIOUS  TIME  STEP 
400  IF  (NITR  .EQ.  0)  THEN 
GAMK(l)  »  GAMMA(l) 

GAMK(2)  =  GAMMAC2) 

ENDIF 

DO  125  L  -  l.NAIRFO 

LI  *  (L-1)*N0DT0T 

KI  -  (L-1)*NP1 

KK  -  (L-l)*2 

DO  130  K  ■  1,2 

IF  (K  .EQ.  1)  I  «  1 

IF  (K  .EQ.  2)  I  *  NODTOT 

XMID  »  0.5  *  CXICI+KI)  +  XICI+KI+I)) 

YMID  -  0.5  *  (YKI+KI)  +  YKI+KI+l)) 

VTANG  =  ( ( 1.  +UG( I+LI ) )*C0SALF( L) -VG( I+LI )*SINALF( L)+UX( L) ) 
+*COSTHL( I+KI )+  ( ( 1.  +UG( I+LI ) )*SINALF(L)+VG( I+LI)*COSALF(L)+UY( L) ) 
+*SINTHL( I+KI )  +  OM£GA( L)*( YMID*COSTHL( I+KI ) 

+-  XMID*SINTHL( I+KI ) ) 

AA(K+KK, 1)  *  -  AAN( I+LI ,NP3)*SS( 1)/WAKE( 1) 

AA( K+KK , 2 )  *  -  AAN(I+LI ,NP4)*SS( 2)/WAKE( 2) 

BB(K+KK)  *  VTANG  +  AAN( I+LI ,NP3)*SS( 1)*GAMMA( 1)/WAKE( 1)  + 

+  AAN(I+LI,NP4)*SS(2)*GAMMA(2)/WAKE(2) 

DO  119  J  «  1, NODTOT 

AA( K+KK , 1 )  *  AA( K+KK , 1 )  +  AAN(I+LI,J)  -  BBN( I+LI , J)*B1( J) 

AA(K+KK,2)  «  AACK+KK.2)  +  AAN(I+LI,J+NODTOT)  -  BBN( I+LI , J)*B2( J) 
BB(K+KK)  »  BBC K+KK)  -  BBN( I+LI , J)*B3( J) 

119  CONTINUE 

DO  120  JJ  «  NODTOT+1 , 2+NODTOT 

AACK+KK.l)  =  AA(K+KK, 1)  -  BBN( I+LI , JJ)*B 1 ( JJ) 

AA(K+KK,2)  *  AACK+KK.2)  -  BBNCI+LI , JJ)*B2C JJ) 

BBC K+KK)  ■  BBC K+KK)  -  BBNC I+LI , JJ)*B3C JJ) 

120  CONTINUE 

ADD  CORE  VORTEX  CONTRIBUTION 

IF  CM.  LE.  1)  GOTO  100 

BBCK+KK)  =  BBC K+KK)  +  SUMCCTCI+LI) 

100  CONTINUE 
130  CONTINUE 
125  CONTINUE 

IF  CNGIES  .EQ.  1)  GOTO  145 

SATISFYING  KUTTA  CONDITION  —  SOLVE  FOR  VORTEX  STRENGTH 

DO  135  I  ■  1,2 
LI  «  CI-1)*NAIRF0+1 
AAACI)  «  AAC LI , 1 )**2-AAC LI+1 , 1 )**2 
BBBCI)  «  AA(LI,2)**2-AACLI+1,2)**2 

CCCCI)  «  2*CAACLI,1)*BBCLI)-AACLI+1,1)*BBCLI+1)-C2-I)*SSC1)/DT) 
DDDCI)  »  2*CAACLI,2)*BBCLI)-AACLI+1,2)*BBCLI+1)-CI-1)*SSC2)/DT) 
EEECI)  »  2*C  AAC  LI , 1 )*AAC  LI , 2 ) -AAC  LI+1 , 1 )*AAC  LI+1 , 2 ) ) 

FFFCI)  *  BBCLI)**2-BBCLI+1)**2+2*SSCI)*GAMMACI)/DT 
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135  CONTINUE 
60  DO  140  I  *  1,2 

DGAMK(I,1)  *  2*AAA( I)*GAMK( 1)+EEE( I)*GAMK( 2)+CCC( I ) 

DGAMK(I,2)  =  2*BBB(I)*GAMK(2)+EEE(I)*GAMK(1)+DDD(I) 

140  DGAMK(I,3)  *  -(AAA(I)*GAMK(1)**2+CCC(I)*GAMK(1)+BBB(I)* 

+  GAMK( 2 )**2+DDD( I )*GAMK( 2 )+EEE ( I ) *GAMK( 1 ) *GAMK( 2 ) +FFF( I ) ) 

C140  WRITE  (6,*) 'DGAMK' , I ,DGAMK( 1,1) ,DGAMK( I ,2) ,DGAMK( I ,3) 

DO  144  1-1,2 
DO  144  KKK=1 , 3 

144  COEFL(I.KKK)  -  DGAMK ( I ,KKK) 

C 

C  GAUSSIAN  ELIMINATION 
C 

R=DGAMK( 2 , 1) /DGAMK ( 1,1) 

DO  200  K  -  2,3 

200  DGAMK(2,K)  *  DGAMK ( 2 ,K) *R*DGAMK( 1 ,K) 

DG(2)  =  DGAMK(2,3)/DGAMK(2,2) 

DG( 1)  *  (DGAMK( 1 , 3) -DGAMK( 1, 2)*DG( 2) )/DGAMK( 1 , 1) 

GAMK(l)  *  GAMK( 1)+DG( 1) 

GAMK(2)  =  GAMK(2)+DG( 2) 

IF  ((ABS(DG( 1)/GAMK( 1))  .  LT. .  0001)  .AND. 

+  ( ABS(DG( 2)/GAMK( 2) ).  LT. .  0001))  GO  TO  300 

COUNT  *  COUNT  +  1 
IF  (COUNT  .EQ.  50)  GO  TO  310 
GO  TO  60 

310  WRITE(6,*)'KUTTA  CONDITION  NOT  SATISFIED  —  COUNT  EXCEEDED  50' 
RETURN 
C 

C  SET  UP  KUTTA  CONDITIONS  FOR  GIESING  CONDITION 
C 

145  DO  451  1*1,2 

LI  *  ( I-1)*NAIRF0+1 

DGAMK(I.l)  *  AA(LI , 1)  +  AA(LI+1,1) 

DGAMK ( 1,2)  *  AA( LI ,2)  +  AA(LI+1,2) 

451  DGAMK( 1,3)  -  -(BB(LI)  +  BB(LI+1)) 

C 

C  SOLVE  KUTTA  CONDITION  BY  GAUSS 
C 

R=DGAMK( 2,1) /DGAMK( 1,1) 

DO  461  K  *  2,3 

461  DGAMK(2,K)  =  DGAMK(2,K)-R*DGAMK(1,K) 

GAMK( 2)  *  DGAMK(2,3)/DGAMK(2,2) 

GAMK(l)  *  (DGAMK( 1 , 3) -DGAMK( 1 , 2)*GAMK( 2) ) /DGAMKC 1,1) 

C 

C  CALCULATE  SOURCE  STRENGTH 
C 

300  DO  160  L  -  l.NAIRFO 
LI  =(L-l)*NODTOT 
DO  160  I  *  l.NODTOT 

160  QK(I+LI)  *  GAMK( 1)*B1(I+LI)  ♦  GAMK(2)*B2(I+LI)  +  83(I+LI) 

C 

C  CALCULATE  TANGENTIAL  VELOCITY  AT  THE  TRAILING  EDGE  BY  BACK 
C  SUBSTITUTION.  THE  PRODUCT  OF  THE  TRAILING  EDGE  VELOCITY 
C  SHOULD  BE  NEGATIVE  .  THIS  IS  CHECKED  IN  THE  MAIN  PROGRAM. 

C 

VTAG1  *  AA(1,1)*GAMK(1)+AA(1,2)*GAMK(2)+BB(1) 
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VTAG50  *  AA( 2 , 1 )*GAMK( 1 )+AA( 2 , 2 )*GAMK( 2 )+BB( 2 ) 

VTAG51  =  AA( 3 , 1)*GAMK( 1)+AA( 3 , 2)*GAMK( 2)+BB( 3) 

VTA100  »  AA(4,1)*GAMK(1)+AA(4,2)*GAMK(2)+BB(4) 

PVTAG  «  VTAG1*VTAG50 
2222  FORMAT( IX, * VTAG* ,4F10. 6) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  SUBROUTINE  TEWAK  (SINALF.COSALF) 

C 

C  COMPUTE  WAKE  ELEMENT  AT  THE  TRAILING  EDGE 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  VYVK(L)  _  Y-VELO  COMP  AT  THE  MID-POINT  OF  THE  WAKE  ELEMENT 

C  WITH  RESPECT  TO  THE  LOCAL  FROZEN  FRAME  OF  REFERENCE 

C  VXWK(L)  _  X-VELO  COMP  AT  THE  MID-POINT  OF  THE  WAKE  ELEMENT 

C  WITH  RESPECT  TO  THE  LOCAL  FROZEN  FRAME  OF  REFERENCE 

C  _ 

SUBROUTINE  TEWAK 

COMMON  /BOD/  IFLAG,NL0WER,NUPPER,N0DT0T,X(202) ,Y(202) , 

+  COSTHEC201) ,SINTKE(201) ,SS(2) ,NP1 ,NP2,NP3,NP4, 

+  NP5 ,XSHIFT, YSHIFT,NAIRFO,XI( 202) ,YI(202) , 

+  COSTHL(201) ,SINTHL(201) 

COMMON  /COF/  A( 201, 211) ,KEQNS 

COMMON  /SING/  Q( 200) ,GAMMA( 2) ,QK(200) ,GAMK(2) 

COMMON  /WAK/  VYW(2) ,VXW(2) ,WAKE(2) ,DT 
COMMON  /WAK2/  VYWK(2) ,VXWK(2) 

COMMON  /CORV/  CV(2 ,200)  ,XC( 2,200)  ,YC(2, 200)  ,M,TD,CWX(2, 200)  , 

+  CWY(2,200)  ,XCI(2,200),YCI(2,200) 

COMMON  /INF1/  AAN( 201,201) ,BBN(201,201) ,AYNP1( 2 ,201) ,BYNP1(2 , 201) , 
+  SUMAAN( 201 , 2) ,SUMBBN( 201 ,2) 

COMMON  /INF2/  SUMCCN( 201 ) ,SUMCCT( 201) ,CYNP1( 2,200) , 

+  CXNP1(2,200) 

COMMON  /GEOM/  SINALF(2),C0SALF(2),0MEGA(2),UX(2),UY(2), PIVOT, 

+  XPRM.YPRM 

COMMON  /GUST/  UG(200) ,VG( 200) ,XGF,UGUST,VGUST 
DIMENSION  C0NT1X( 2 ) , CONTI Y( 2 ) , CONT2XC  2 ) , CONT2Y( 2 ) , 

+  CONT3XC  2 ) , CONT3YC  2 ) , C0NT4X( 2 ) , C0NT4Y( 2 ) , CONT5X( 2 ) , CONT5Y( 2 ) 

DO  20  I  -  l.NAIRFO 
NN  *  (I-1)*NP1 

XMID  *  0.5  *  (XCNP1+NN)  +  X(NP2+I-1)) 

YMID  «  0.5  *  (Y(NP1+NN)  +  Y(NP2+I-1)) 

UGW  =  0.  0 

VGW  «  0.0 

XG  =  XMID 

IF  (XG  .GT.  XGF)  GO  TO  10 

UGW  *  UGUST 

VGW  »  VGUST 

10  VYWK(I)  »  ( 1. +UGW)*SINALF( I )+VGW*COSALF( I) 

VXWK(I)  »  ( 1. +UGW)*COSALF( I ) -VGW*SINALF( I ) 

2221  FORMAT( IX, ’SINALF  :  ' ,4F10.  5) 

CONTIY(I)  »  ( 1.  +UGW)*SINALF( I )+VGW*COSALF( I ) 

CONTIX(I)  «  ( 1.  +UGW)*COSALF(I)-VGW*SINALF(I) 

CONT2Y( I)  «  0 
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C0NT2X( I)  =  0 
C0NT3Y( I)  =  0 
C0NT3X( I )  =  0 
C 

C  CONTRIBUTION  FROM  THE  SOURCE  AND  WAKE  DIST  FROM  THE  AIRFOIL 
C 

DO  110  L  *  l.NAIRFO 
KJ  «  (L-1)*N0DT0T 
DO  120  J  ■  l.NODTOT 

VYWK(I)  ■  VYVK(I)  +  AYNP1(I,J+KJ)*QK(J+KJ)  +  BYNP1(I,J4-KJ)*GAMK(L) 
VXWK(I)  »  VXWK(I)  -  BYNPK I , J+KJ)*QK( J+KJ)  +  AYNP1(I, J+KJ)*GAMK(L) 
C0NT2Y(I)  «  C0NT2Y(I)+  AYNP1(I,J+KJ)*QK(J+KJ) 

C0NT2X( I )  *  C0NT2X(I)-  BYNP1( I, J+KJ)*QK( J+KJ) 

C0NT3Y(I)  »  C0NT3Y(I)+  BYNP1C I , J+KJ)*GAMK(L) 

120  C0NT3XCI)  «  C0NT3X(I)+  AYNP1(I, J+KJ)*GAMK(L) 

110  CONTINUE 
C 

C  ADD  WAKE  ELEMENT  CONTRIBUTION  OF  THE  OTHER  AIRFOIL 
C 

J  «  3-1 

VYWK(I)  =  VYVK(I)+  BYNP1(I ,NP3)*(GAMMA( J)-GAMK( J)) 

+  *  SS( J)/WAKE( J) 

VXWK(I)  =  VXWK(I)+  AYNP1(I,NP3)*(GAMMA( J)-GAMK( J)) 

+  *  SS( J)/WAKE( J) 

C0NT4Y( I)  »  BYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 

+  *  SS( J)/WAKE( J) 

125  C0NT4X( I )  »  AYNP1(I,NP3)*(GAMMA(J)-GAMK(J)) 

+  *  SS( J)/WAKE( J) 

C 

C  ADD  CORE  VORTEX  CONTRIBUTION 

C 

C0NT5Y(I)  *  0 

C0NT5X( I)  =  o 
IF  (M  .EQ.  1)  GO  TO  140 
MM1  =  M  -  1 

DO  150  L  *  l.NAIRFO 
KN  =(L-1)*MM1 

DO  130  N  =  1.MM1 

VYWK(I)  =  VYWK(I)  +  CYNP1( I ,N+KN)*CV(L,N) 

VXWK(I)  »  VXWK(I)  +  CXNP1(I ,N+KN)*CV(L,N) 

C0NT5Y(I)  *  CYNP1(I,N+KN)*CVCL,N) 

130  C0NT5X(I)  =  CXNP1( I ,N+KN)*CV(L,N) 

150  CONTINUE 
140  CONTINUE 
20  CONTINUE 

2222  FORMAT  (IX,  ’VXWK  :',6F10.5) 

2223  FORMAT  (IX,  ’VYWK  :',6F10. 5) 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  c 

C  SUBROUTINE  PRESS  (SINALF, COS ALF, OMEGA, UX.UY)  C 

C  C 

C  COMPUTE  UNSTEADY  FLOW  PRESSURE  DISTRIBUTION  C 

C  AND  VELOCITY  POTENTIAL  AT  MID-POINTS  OF  PANELS  C 
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CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


VTANG  .  TANGENTIAL  VELO  COMP  OF  A  FIXED  POINT  OF  THE  MOVING 

LOCAL  FRAME  OF  REFERENCE. 

PHIK(I) .  DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 

I-TH  PANE  AT  THE  CURRENT  TIME  STEP  FOR  THE  UNSTEADY 
FLOV  CASE. 

PHI(I)  .  DISTURBANCE  VELO  POTENTIAL  AT  THE  MID  POINT  OF  THE 

I-TH  PANEL  FOR  THE  PREVIOUS  TIME  STEP 


PHILE(L). . . .  DIFFERENCE  OF  THE  POTENTIALS  OF  THE  LEADING  EDGE  TO 
THE  LOWER  TRAILING  EDGE  FOR  THE  RESPECTIVE  AIRFOIL 
PINK(L)  ....  DIFFERENCE  OF  THE  POTENTIALS  AT  A  POINT  100  CHORD 

LENGTH  UPSTREAM  OF  THE  LE  FOR  THE  RESPECTIVE  AIRFOIL 

SUMC(L)  -  GAMMA  ASSOCIATED  WITH  THE  INTEGRATION  OF  THE  DISTURB 

ANCE  VELOCITY  AROUND  THE  WHOLE  AIRFOIL 

THE  FOLLOWING  INFLUENCE  COEF  ARE  COMPUTED  FOR  A  FINER  GRID  ON  THE 
AIRFOIL  SO  AS  TO  OBTAIN  A  MORE  ACCURATE  VELO  POTENTIAL  AT  THE  TE 
THE  INFLUENCE  COEF  ON  THE  I-TH  PANEL  FROM  THE  J-TH  PANEL  OF  THE 
AIRFOIL  REMAINS  THE  SAME  FOR  ALL  TIME  STEP 


AANP1C I , J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 
OF  AIRFOIL  1  DUE  TO  UNIT  STRENGTH  DIST  SOURCE  STRENG¬ 
TH  ON  THE  J-TH  PANEL  OF  AIRFOIL  1 

BBNP1( I , J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 
OF  AIRFOIL  1  DUE  TO  UNIT  STRENGTH  DIST  VORTICITY  STR¬ 
ENGTH  ON  THE  J-TH  PANEL  OF  AIRFOIL  1 

AANP2( I , J,K).  .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 
OF  AIRFOIL  2  DUE  TO  UNIT  STRENGTH  DIST  SOURCE  STRENG¬ 
TH  ON  THE  J-TH  PANEL  OF  AIRFOIL  2 

BBNP2( I , J,K). .  NORMAL  VELOCITY  INDUCED  AT  THE  I-TH  PANEL  SUB  NODE  K 
OF  AIRFOIL  2  DUE  TO  UNIT  STRENGTH  DIST  VORTICITY  STR¬ 
ENGTH  ON  THE  J-TH  PANEL  OF  AIRFOIL  2 


SUBROUTINE  PRESS 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHE(201) .SINTHEC 201) ,SS( 2) ,NP1 ,NP2 ,NP3,NP4, 

+  NP5 ,XSHIFT,YSHIFT,NAIRF0,XI(202) ,YI(202) , 

+  COSTHL(201) ,SINTHL(201) 

COMMON  /CPD/  CP(200),SCL,T,SCM,SGAM 
COMMON  /NUM/  PI.PI2INV 

COMMON  /SING/  Q(200) ,GAMMA(2) ,QK( 200) ,GAMK(2) 

COMMON  /WAX/  VYW(2) ,VXW(2) ,WAKE(2) ,DT 

COMMON  /CORV/  CV(2,200),XC(2,200),YC(2,200),M,TD,CWX(2,200), 

+  CWY(2,200)  ,XCI( 2,200)  ,YCI(2,200) 

COMMON  /INF1/  AAN( 20 1,201) ,BBN(201,201) ,AYNP1(2,201) ,BYNP1(2,201) , 
+  SUMAAN(201,2) ,SUMBBN(201,2) 

COMMON  /INF2/  SUMCCN(201) ,SUMCCT(201) ,CYNP1(2,200) , 

♦  CXNP1(2,200) 

COMMON  /INFL3/  AANP(101,101,6),BBNP(101,101,6) 

COMMON  /POT /  PHI( 200) ,PHIK( 200) 

COMMON  /GUST/  UG( 200 ) , VG( 200 ) , XGF , UGUST , VGUST 
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COMMON  /EXTV/  UE(200) ,VN(200) 

COMMON  /BMAT/  B1,B2,B3 

COMMON  /GEOM/  SINALF(2) ,C0SALF(2) ,0MEGA(2) ,UX(2) ,UY(2) , PIVOT, 

+  XPRM.YPRM 

DIMENSION  PHITEL( 2 ) , PHITEU( 2 ) ,WGHT( 5 ) , PLOC( 5 ) , SUMC( 2 ) , 
+AANP1(50,50,9) ,AANP2(50,50,9) ,BBNP1(50,50,9) ,BBNP2(50 ,50 ,9) , 
+C0STHPC 102,6) , SINTHPC 102,6), AANP4( 2 ) , PHIL( 102 ) , UGU( 100,6), 

+VGU( 100,6) ,PDUM(2) ,DP(2) ,BBNP4(2) ,B1( 200) ,B2(200) ,B3(200) ,PHILE(2) 
3333  DIMENSION  C0N2(2) ,VTGT(2,50) ,WAKCON(2,50) ,VOTCON(2,50) , 

+  DUMC0N(2,50) 

7733  DIMENSION  CP1(200) 

4444  REAL  *8  PINK(2),VELX 

DATA  WGHT/.  11846344,.  23931434,.  28444444, 

+  .23931434,. 11846344/ 

DATA  PLOC/.  04691008 , .  23076535 , . 50000000 , 

+  .76923466,. 95300899/ 

WRITE  (6,1000) 

IF  (M  .GT.  1)  GO  TO  510 


COMPUTE  THE  INFLUENCE  COEFFICIENT  FOR  A  FINER  GRID 

_  ON  THE  AIRFOIL  BY  THE  SAME  AIRFOIL  . . . 

(  ONLY  COMPUTED  ONCE  ) 


100 


120 

200 

510 


DO  200  L  »  1 ,NAIRFO 
LI  =  (L-l)*NODTOT 

KI  =  (L-1)*NP1 

DO  200  I  =  1 , NODTOT 
DO  200  K  =  1,5 

DX  =  PLOC( K)*( X( I+KI+1 ) -X( I+KI ) ) 

DY  =  PLOC(K)*(Y( I+KI+1 )-Y( I+KI)) 

XMID  =  X( I+KI )  +  DX 

YMID  *  Y( I+KI)  +  DY 

DO  200  J  =  1, NODTOT 

FLOG  =0.0 
FT AN  =  PI 

IF  ( J+LI  .EQ.  I+LI)  GO  TO  100 
DXJ  =  XMID  -  X( J+KI) 

DXJP  =  XMID  -  X( J+KI+1) 

DYJ  =  YMID  -  Y( J+KI) 

DYJP  =  YMID  -  Y( J+KI+1) 

FLOG  =  .  5*ALOG((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
FTAN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE( I+KI )*C0STHE( J+KI)  + 

+  SINTHEC I+KI )*SINTHE( J+KI) 

STIMTJ  =  SINTHEC I+KI )*COSTHE( J+KI)  - 
+  COSTHEC I+KI )*SINTHE( J+KI) 

IF  (L  .EQ.  2)  GOTO  120 

AANP1(I,J,K)  =  PI2INV*( FT AN* CTIMTJ  +  FLOG+STIMTJ) 

BBNP1( I , J,K)  =  PI2INV*(FLOG*CTIMTJ  -  FTAN*STIMTJ) 

GOTO  200 

AANP2(I,J,K)  «  PI2INV*(FTAN*CTIMTJ  +  FLOG*STIMTJ) 
BBNP2(I,J,K)  -  PI2INV*( FLOG* CTIMTJ  -  FTAN+STIMTJ) 

CONTINUE 

CONTINUE 


C 
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FIND  TANGENTIAL  VELOCITY  AT  THE  MIDPOINT  OF  THE  I-TH  PANEL 

DO  600  L  =  1.NAIRF0 
LI*( L- 1 )*N0DT0T 
KI=(L-1)*NP1 
SUMC(L)  =0.0 
DO  600  I  =  1 , NODTOT 
C0NTR1  =0.0 

C0NTR2  =0.0 

VAKC0N(L,I)  =0.0 
VOTCON(L.I)  =  0.0 
DUMCON(L.I)  =  0.0 

XMID  =  0.5  *  (XICI+KI)  +  XI( I+KI+1)) 

YMID  =  0.  5  *  (YI(I+KI)  +  Ylfl+KI+D) 

DX  =  ( XI  ( I+KI+1)  -  XKI+KI)) 

DY  =  (YI( I+KI+1)  -  YI( I+KI) ) 

DIST  =  SQRT( DX*DX+DY*DY ) 

ACCOUNT  FOR  THE  FREESTREAM  AND  THE  MOTION  OF  THE  AIRFOIL 

VSX  =  ( 1.  +UG( I+LI) )+C0SALF(L) -VG( I+LI)*SINALF(L)  +  OMEGA(L)* 

+  YMID  +  UX( L) 

VSY  =  (l.+UG(I+LI))*SINALF(L)+VG(I+LI)*COSALF(L) 

+  -  OMEGA(L)*XMID  +  UY(L) 

VS  =  VSX*VSX  +  VSY*VSY 

VNORM  =  -VSX*SINTHL( I+KI )+VSY*COSTHL( I+KI) 

VTANG  =  ( ( ( 1.  +UG( I+LI) )*COSALF(L)-VG( I+LI)*SINALF(L)+UX(L)) 

+  *COSTHL( I+KI )+  ((1.  +UG(I+LI))*SINALF(L)+VG(I+LI)*COSALF(L)+UY(L)) 
+  *SINTHL( I+KI)+  OMEGA( L)*( YMID*COSTHL( I+KI ) 

+  -  XMID*SINTHL( I+KI ) ) ) 

VTFREE  =  VTANG 

VACT  =  VTANG 

INTRODUCE  SMALLER  GRIDS  FOR  THE  PURPOSE  OF  THE  VELO  POTENTIAL. 

VELO  ONLY  CALCULATED  AT  THE  MI DPT  OF  THE  PANEL  WHERE  K  =  3 

DO  260  K  =  1,5 

DX  =  PLOC(K)*(X( I+KI+1 )-X( I+KI)) 

DY  =  PLOC(K)*(Y( I+KI+1)-Y( I+KI)) 

XINT  =  X( I+KI )  +  DX 
YINT  =  Y( I+KI)  +  DY 
VDUM  =0.0 

INFLUENCE  COEF  ON  I-TH  PANEL  DUE  TO  THE  WAKE  ELEMENTS 

DO  245  MM  =  1,  NAIRFO 
J  =  NP1+MM 

DXJ  =  XINT  -  X(J) 

DXJP  =  XINT  -  X(NAIRF0*NP1+MM) 

DYJ  =  YINT  -  Y(J) 

DYJP  =  YINT  -  Y(NAIRF0*NP1+MM) 

FLOG  «  . 5*ALOG( ( DXJP*DXJP+DYJP*DYJP) / ( DXJ*DXJ+DYJ*DYJ) ) 

FT AN  =  ATAN2( DYJP+DXJ -DXJP+DYJ ,DXJP*DXJ+DYJP*DYJ) 

CTIMTJ  =  COSTHE( I+KI)*COSTHE( J)  + 

+  SINTHE(I+KI)*SINTHE( J) 

STIMTJ  =  SINTHEC I+KI)*COSTHE( J)  - 
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+  COSTHE( I+KI)*SINTHE( J) 

AANP4(MM)  =  PI2INV*(FTAN*CTIMTJ  +  FLOG+STIMTJ) 

245  BBNP4(MM)  =  PI2INV*(FL0G*CTIMTJ  -  FTAN*STIMTJ) 

C0NTR1  =  SS(1)*(GAMMA(1)-GAMK(1))*AANP4(1) 
+/WAK£(1)+SS(2)*(GAMMA(2)-GAMK(2))*AANP4(2)/WAKE(2) 

CONTRIBUTION  TO  VELO  COMPONENT  BY  WAKE  ELEMENT. 

IF  (K  .EQ.  3)  THEN 

VACT  =  VACT  +  SS(1)*(GAMMA(1)-GAMK(1))*AAN(I+LI, 
+NP3)/WAKE(1)+SS(2)*(GAMMA(2)-GAMK(2))*AAN(I+LI,NP4)/WAKE(2) 
VNORM  *  VNORM  +  SS( 1)*(GAMMA( 1)-GAMK( 1) )*BBNP4( 1) 

+/WAKE( 1)+SS( 2)*(GAMMA( 2) -GAMK( 2) )*BBNP4( 2)/WAKE( 2) 

ENDIF 


EFFECTS  ON  AIRFOIL  BY  THE  SAME  AIRFOIL 


INTEGRATION  AROUND  FIRST  AIRFOIL 

CONTRIBUTION  TO  VELO  COMP  BY  AIRFOIL  1  WHEN  K  =  3 

DO  300  J  =  l.NODTOT 
IF  (L  .EQ.  2)  GOTO  270 

VDUM  =  VDUM-BBNP1( I , J,K)*QK( J)+AANP1( I , J,K)*GAMK( 1) 

IF  (K  .EQ.  3)  THEN 

VACT  -  VACT-BBN( I+LI , J)*QK( J)+ 

+AAN( I+LI , J)*GAMK( 1) 

VNORM  =  VN0RM+(AANP1( I , J,K)*QK( J))+ 

+(BBNP1( I , J,K)*GAMK( 1) ) 

ENDIF 
GOTO  300 

INTEGRATION  AROUND  SECOND  AIRFOIL 

CONTRIBUTION  TO  VELO  COMP  BY  AIRFOIL  2  WHEN  K  -  3 

270  VDUM  =  VDUM-BBNP2C I , J,K)*QK( J+N0DT0T)+AANP2( I , J,K)*GAMK(2) 
IF  (K  .EQ.  3)  THEN 

VACT  =  VACT-BBN( I+LI , J+NODTOT)*QK( J+NODTOT)+ 

+  AAN( I+LI , J+NODTOT ) *GAMK( 2 ) 

VNORM  =  VN0RM+(AANP2( I ,J,K)*QK( J+NODTOT) )+ 

+  ( BBNP2( I , J,K)*GAMK( 2) ) 

ENDIF 

300  CONTINUE 


EFFECTS  ON  AIRFOIL  BY  THE  OTHER  AIRFOIL 


INTEGRATION  AROUND  BOTH  AIRFOIL 

CONTRIBUTION  TO  VELO  COMP  BY  BOTH  AIRFOIL  WHEN  K  ■  3 

DO  350  J  ■  NODTOT+2 , 2+NODTOT+l 
DXJ  «  XINT  -  X(J-KI) 

DXJP  ■  XINT  -  X( J-KI+1) 

DYJ  =  YINT  -  Y(J-KI) 

DYJP  =  YINT  -  Y( J-KI+1) 

FLOG  *  .5*AL0G((DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
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nan  on 


FTAN  =  ATAN2  C D Y JP*DX J -DX JP*DY J , DX JP*DX J+DY JP*DY J ) 
CTIMTJ  =  COSTHE( I+KI )*COSTHE( J-KI )  + 

+  SINTHEC I+KI)*SINTHE( J-KI) 

STIMTJ  =  SINTHEC I+KI )*COSTHE( J-KI)  - 
+  COSTHEC I+KI )*SINTHEC J-KI) 

AANP3  =  PI2INV*(FTAN*CTIMTJ  +  FLOG+STIMTJ) 

BBNP3  =  PI2INV*(FLOG*CTIMTJ  -  FTAN+STIMTJ) 

IF  (K.EQ.  3)  THEN 

VACT  *  VACT-BBN( I+LI , J-LI - 1 )*QK( J-LI - 1 )+ 

+  AAN( I+LI, J-LI- 1)*GAMK(3-L) 

VNORM  *  VNORM+C  AANP3*QK( J-LI - 1 ) )+ 

+  (BBNP3+GAMKC3-L)) 

ENDIF 

350  VDUM  *  VDUM  -  BBNP3*QK( J-LI-1)  +  AANP3*GAMK( 3-L) 

CONTRIBUTION  BY  CORE  VORTICES 

1.  TO  VELO  POTENTIAL 

2.  TO  VELO  COMPONENT  ONLY  WHEN  K  *  3 

IF  (M. EQ. 1)  GOTO  150 
MM1  =  M  -  1 
SUMCN  =  0.  0 
SUMCT  =0.0 
DO  400  MM  =  1,  NAIRFO 
KN  =  (MM-1)*MM1 
DO  400  N  =  1 ,MM1 
DXC  =  XI NT  -  XC(MM.N) 

DYC  *  YINT  -  YC(MM,N) 

DISTC  =  SQRT(DXC*DXC+DYC*DYC) 

COSTHN  =  DXC/DISTC 

SINTHN  =  DYC /DISTC 

CTIMTN  =  COSTHEC I+KI )*COSTHN  +  SINTHEC I+KI)*SINTHN 

STIMTN  =  SINTHEC I+KI )*C0STHN  -  COSTHEC I+KI )*SINTHN 

CCN  =  -CTIMTN/DISTC 
CCT  =  -STIMTN/DISTC 
SUMCN  =  SUMCN  +  CCN*CVCMM,N) 

400  SUMCT  =  SUMCT  +  CCT*CVCMM,N) 

IF  CK  .EQ.  3)  THEN 

VACT  =  VACT+SUMCT*PI2INV 

VNORM  *  VNORM+C PI2INV*SUMCN) 

ENDIF 

CONTR2  *  PI2INV*SUMCT 
150  CONTINUE 

VTANG  ■  VTANG  +  DBLE( CCONTRl+CONTR2+VDUM)*WGHT(K) ) 
SUMCCL)  «  SUMCC L)+VDUM*DIST*WGHTC K) 

WAKCONCL.I)  «  WAKCONCL, I)+C0NTR1 
VOTCONCL.I)  ■  V0TC0NCL,I)+C0NTR2 
DUMCONCL, I)  ■  DUMCON ( L , I ) +VDUM 
260  CONTINUE 

PHIKCI+LI)  «  (( VTANG) -VTFREE)*DIST 
CPC I+LI)  »  VS  -  CVACT*VACT) 

7755  CP1CI+LI)  «  CPC I+LI) 

UECI+LI)  «  VACT 

VNCI+LI)  »  VNORM 

600  CONTINUE 

6688  FORMATC2I5 ,7E15.  7) 


on  n  n  n  o o  n  o  o  o  o  o  o 


3335  FORMAT  (2I5.7F10.  6) 

COMPUTE  DISTURBANCE  POTENTIAL  AT  THE  LEADING  EDGE  BY  LINE 
INTEGRAL  OF  THE  VELOCITY  FIELD 
FROM  UPSTREAM  (AT  INFINITY)  TO  THE  LEADING  EDGE 

TRY  =1 
PINK(l)  =0.0 
PINK(2)  =0.0 
130  DO  56  L  =  l.NAIRFO 

YMID  =  PIVOT  *  SINALF(L)  +  (L-1)*YSHIFT 

XLE  =  -PIV0T*C0SALF(L)+(L-1)*XSHIFT 

XL  =  XLE 

XL  =0.0 

NPHI  =  10  *  NLOWER 

PDUM(L)  =  0.0 

DO  30  I  =  l.NPHI 

FRACT  =  FLO AT ( I ) /FLO AT ( NPH I ) 

XLP  =  -100.0  *TRY  *  (1.0  -  C0S(0.  5*PI*FRACT))+XLE 

IF  (I  .EQ.  1)  XLP  =  -0.  000197+XLE 

XLP  =  -10.0  *  (1.0  -  C0S(  0.  5*PI*FRACT) ) 

DELX  =  XL  -  XLP 

XMID  =  0. 5*(XL+XLP) 

XMID  =  0.  5*(XL+XLP)*C0SALF(L) 

YMID  =  0.5*(XL+XLP)*SINALF(L) 

XL  =  XLP 

VELX  =  UGUST 


ADD  CONTRIBUTION  OF  J-TH  PANEL 


24 


DO  40  MM  =  l.NAIRFO 
LJ  =( MM - 1 )*NODTOT 

KJ  =(MM-1)*NP1 

DO  20  J  =  1,NP1 
IF  (J  .EQ.  NP1)  GO  TO  24 
DXJ  *  XMID  -  X( J+KJ) 

DXJP  =  XMID  -  X( J+KJ+1) 

DYJ  =  YMID  -  Y( J+KJ) 

DYJP  =  YMID  -  Y( J+KJ+1) 

FLOG  1  =  . 5*ALOG( ( DX JP*DXJP+DYJP*DY JP) / ( DXJ*DXJ+DYJ*DYJ) ) 

FT AN  =  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CALMTJ  =  -COSTHE( J+KJ) 

SALMTJ  =  SINTHE( J+KJ) 

CALMTJ  =  -COS ALF ( L) +COSTHE ( J+KJ )  -  SINALF(  L)*SINTHE( J+KJ) 
SAIUTJ  =  -SINALF(L)*COSTHE( J+KJ)  +  COSALF(L)*SINTHE( J+KJ) 
APY  =  PI2INV*( FTAN+CALMTJ  +  FLOG* SALMTJ) 

BPY  =  PI2INV*(FLOG*CALMTJ  -  FTAN*SALMTJ) 

VELX  »  VELX  -  DPROD( BPY , QK( J+LJ) )  +  DPROD( GAMK(MM) , APY) 
GO  TO  20 

DXJ  »  XMID  -  X( J+KJ) 

DXJP  =  XMID  -  X(2*NP1+MM) 

DYJ  =  YMID  -  Y( J+KJ) 

DYJP  =  YMID  -  Y( 2+NP1+MM) 

FLOG  =  .  5*ALOG( ( DXJP*DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 
FTAN  *  ATAN2(DYJP*DXJ-DXJP*DYJ,DXJP*DXJ+DYJP*DYJ) 

CALMTJ  =  -COSALF(L)*COSTHE( J+KJ)  -  SINALF( L)*SINTHE( J+KJ) 
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SALMTJ  «  -SINALF( L)*COSTHE( J+KJ)  +  COSALF(L)*SINTHE( J+KJ) 

APY  =  P I 2 I NV* ( FTAN*C ALMT J  +  FLOG* SALMTJ) 

DUMMY  =  SS(MM)*(GAMMA(MM)-GAMK(MM))/WAKE(MM) 

VELX  *  VELX  +  DPROD( APY, DUMMY) 

20  CONTINUE 

40  CONTINUE 

ADD  CORE  VORTEX  CONTRIBUTION 

IF  (M  .EQ.  1)  GO  TO  50 
MM1  *  M  -  1 
DO  70  II  -  l.NAIRFO 
KN  «  (II-1)*MM1 
DO  60  N  »  1.MM1 
DX  «  XMID  -  XC(II,N) 

DY  «  YMID  -  YC(II,N) 

DIST  *  SQRT( DX*DX+DY*DY) 

COSTHN  =  DX/DIST 
SINTHN  *  DY/DIST 

SALMTN  «  -SINALF( L)*COSTHN  +  COSALF(L)*SINTHN 

CPT  *  -PI2INV*SALMTN/DIST 

60  VELX  »  VELX  +  DPROD(CPT,CV(II,N)) 

70  CONTINUE 

50  CONTINUE 

PDUM(L)  *  PDUM(L)  +  VELX  *  DBLE(DELX) 

7771  FORMAT  (2I5.6E14. 6.2F10. 6) 

30  CONTINUE 

DP(L)  =  PDUM(L)-PINK(L) 

56  CONTINUE 

PINK(l)  *  PDUM(l) 

PINK(2)  *  PDUM( 2 ) 

COMPUTATION  OF  THE  VELOCITY  POTENTIAL  FOR  MIDPOINT  OF  EACH  PANEL 

DO  240  L  *  l.NAIRFO 

LI  *  (L-1)*NODTOT 

PHP  «  -PINK(L) 

BEGIN  WITH  LOWER  SURFACE 

DO  230  I  *  NLOWER.1,-1 
PHC  *  PHP-PHIK( I+LI) 

PHIK(I+LI)  »  0. 5*(PHP+PHC) 

230  PHP  «  PHC 

PHITEL(L)  «  PHC 

RESET  FOR  UPPER  SURFACE 

PHP  -  -PINK(L) 

DO  250  I  »  NLOWER+1 , NODTOT 

PHC  ■  PHP+PHIK( I+LI) 

PHIK(I+LI)  *  0.  5*(PHP+PHC) 

250  PHP  ■  PHC 

PHITEU(L)  «  PHC 
240  CONTINUE 


c 

c 


290 

295 

235 


COMPUTE  CP  AT  MID  POINT  OF  I-TH  PANEL 

DO  295  L«  l.NAIRFO 
LI  =  (L-1)*N0DT0T 

XI  =  (L-1)*NP1 

SUMC(L)  =  SUMC(L)/SS(L) 

DO  290  I  «  l.NODTOT 

XIMID  =  .  5*(XI(I+KI)  +  XI(I+KI+1)) 

XMID  *  . 5*(X(I+KI)  +  X(I+KI+1)) 

YMID  ■  .  5*(Y(I+KI)  +  Y(I+KI+1)) 

CP(I+LI)»  CP(I+LI)  -  2. *( PHIK( I+LI ) -PHI ( I+LI ) ) /DT 
WRITE  (6,1050)  I+LI, XIMID, XMID, YMID, QK( I+LI), GAMK(L) ,CP( I+LI) , 
+  UE( I+LI ),VN( I+LI ),PHIK( I+LI), PHI ( I+LI ),SUMC(L) 

CONTINUE 

WRITE  (6,235)  PINK(L) 

CONTINUE 

FORMAT  (IX/ PHI  AT  LEADING  EDGE  *' .F10.6./) 

rnnuAT/  /  /.v  »  t*  /. v  •vt/tx*  r v  *v/r\t  zv  fv/’rx*  /v  1 Q(J)f  6X 


SUBROUTINE  CORVOR  ( S INALF , COS ALF ) 

COMPUTE  THE  LOCAL  VELOCITIES  OF  CORE  VORTICES 


1000  F0RMAT(/,4X, ' J’ ,4X, ’XI( J) ' ,5X, ’X(  J) ’ ,6X, *Y(  J) ’ ,6X, ’Q(* 

+' GAMMA’ ,4X. 'CP( J) 1 ,7X, ' V( J) * ,6X, ' VN( J) ' ,5X, 'PHIK' ,6X, 

+’PHI' ,  6X, ' INTGAMKf ,/) 

1050  F0RMAT( 15 , 12F10.  5) 

1200  F0RMAT( IX, 'LENGTH  OF  LEADING  EDGE  INTEGRATION  IN  CHORD  *’,F10.6/) 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 
C 
C 
C 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  CCVX(I.N)  ....  X-VELO  OF  THE  I-TH  AIRFOIL,  N-TH  CORE  VORTEX  WITH 
C  RESPECT  TO  THE  CURRENT  FROZEN  FRAME  OF  REFERENCE 

C  CCVY(I,N)  _  Y-VELO  OF  THE  I-TH  AIRFOIL,  N-TH  CORE  VORTEX  WITH 

C  RESPECT  TO  THE  CURRENT  FROZEN  FRAME  OF  REFERENCE 

C  XC( I ,N)  _  X-COORD  OF  THE  LOCATION  OF  THE  I-TH  AIRFOIL,  N-TH 

C  CORE  VORTEX  W.  R.  T.  THE  GLOBAL  FRAME  OF  REFERENCE. 

C  YC(I,N)  _  Y-COORD  OF  THE  LOCATION  OF  THE  I-TH  AIRFOIL,  N-TH 

C  CORE  VORTEX  W.  R.  T.  THE  GLOBAL  FRAME  OF  REFERENCE. 

C  GAMY(I)  _  GLOBAL  Y-VELO  AT  A 

C  CORE  VORTEX  DUE  TO 

C  ON  ONE  PANEL 

C  GBMY(I)  _  GLOBAL  Y-VELO  AT  A 

C  CORE  VORTEX  DUE  TO 

C  ON  ONE  PANEL 

C  AMY(I)  _  LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRENGTH 

C  ON  ONE  PANEL 

C  BMY(I)  _  LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  VORTEX  DIST  OF  UNIT  STRENGTH 

C  ON  ONE  PANEL 

C  SUMAMY(I)  _  LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 

C  CORE  VORTEX  DUE  TO  A  SOURCE  DIST  OF  UNIT  STRENGTH 

C  ON  ALL  PANELS 


LOCATION  OF  THE  I-TH  AIRFOIL 
A  SOURCE  DIST  OF  UNIT  STRENGTH 

LOCATION  OF  THE  I-TH  AIRFOIL 
A  VORTEX  DIST  OF  UNIT  STRENGTH 
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n  n  o  o  o 


ooo  non  oooonononnooo 


SUMBMY(I)  .. 


LOCAL  Y-VELO  AT  A  LOCATION  OF  THE  I-TH  AIRFOIL 
CORE  VORTEX  DUE  TO  A  VORTEX  DIST  OF  UNIT  STRENGTH 
ON  ALL  PANELS 

GLOBAL  X-VELO  AT  A  LOCATION  OF  A  CORE  VORTEX  DUE 
TO  ANOTHER  POINT  VORTEX  OF  UNIT  STRENGTH 
GLOBAL  Y-VELO  AT  A  LOCATION  OF  A  CORE  VORTEX  DUE 
TO  ANOTHER  POINT  VORTEX  OF  UNIT  STRENGTH 
LOCAL  X-VELO  AT  A  LOCATION  OF  A  CORE  VORTEX  DUE 
TO  ANOTHER  POINT  VORTEX  OF  UNIT  STRENGTH 
LOCAL  Y-VELO  AT  A  LOCATION  OF  A  CORE  VORTEX  DUE 
TO  ANOTHER  POINT  VORTEX  OF  UNIT  STRENGTH 


SUBROUTINE  CORVOR 

COMMON  /BOD/  IFLAG,NLOWER,NUPPER,NODTOT,X(202) ,Y(202) , 

+  COSTHE(201),SINTH£(201),SS(2),NP1,NP2>NP3,NP4, 

+  NP5 , XSHIFT, YSHIFT , NAIRFO , XI ( 202) , YI ( 202 ) , 

+  COSTHL( 201) ,SINTHL( 201) 

COMMON  /SING/  Q(200) ,GAMMA(2) ,QK( 200) ,GAMK(2) 

COMMON  /WAK/  VYW(2) ,VXW(2) ,WAKE(2) ,DT 

COMMON  /CORV/  CV(2,200) ,XC( 2,200) ,YC(2, 200) ,M,TD,CCVX(2,200) , 
+  CCVY(2,200),XCI(2,200),YCI(2,200) 

COMMON  /POT/  PHI( 200) ,PHIK( 200) 

COMMON  /GUST/  UG(200) ,VG(200) ,XGF,UGUST,VGUST 
COMMON  /NUM/  PI.PI2INV 

COMMON  /GEOM/  SINALF(2) ,COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2) .PIVOT, 
+  XPRM.YPRM 

DIMENSION  AMY( 2) ,BMY(2) , SUMAMY( 2 ) , SUMBMY( 2 ) , 

+VX( 2 ) , VY( 2 ) , GAMYC 2 ) , GBMY( 2 ) 

IF  (M.EQ.  1)  GOTO  40 
MM1  =  M  -  1 

ACCOUNT  FOR  THE  FREE  STREAM  INCLUDING  GUST  EFFECTS 

UGC  *  0.  0 
VGC  *  0.0 
DO  15  I  *  1, NAIRFO 
KN  =  ( I-1)*MM1 

DO  10'  N  -  1.MM1 

XG  -  XC(I,N)*COSALF(I)  +  YC(I,N)*SINALF(I) 

IF  (XG  . GT.  XGF)  GO  TO  5 
UGC  «  UGUST 
VGC  »  VGUST 
CONTINUE 

VY( I)  »  (l.+UGC)*SINALF(I)+VGC*COSALF(I) 

VX(I)  *  ( 1.  +UGC)*COSALF(I)-VGC*SINALF( I) 

XMID  *  XC(I,N) 

YMID  «  YC(I,N) 

CALCULATE  THE  INFLUENCE  COEFFICIENT  DUE  TO  THE  AIRFOILS 

DO  25  L  »  1, NAIRFO 
LJ  «(L-1)*NODTOT 

KJ  =( L-1)*NP1 

SUMAMY(I)  =0.0 
SUMBMY(I)  =  0.  0 
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DO  20  J  «  1 ,NP1 
DXJ  ■  XMID  -  X( J+KJ) 

DYJ  *  YMID  -  Y( J+KJ) 

IF  (J  .EQ.  NP1)  00  TO  11 
DXJP  =  XMID  -  X( J+KJ+1) 

DYJP  *  YMID  -  Y( J+KJ+1) 

GO  TO  12 

11  DXJP  »  XMID  -  X(2*NP1+L) 

DYJP  -  YMID  -  Y(2*NP1+L) 

12  FLOG  »  .  5*AL0G( (DXJP*DXJP+DYJP*DYJP)/(DXJ*DXJ+DYJ*DYJ)) 
FTAN  -  ATAN2( DY JP*DXJ -DXJP*DY J , DXJP+DX J+DY JP*DY J ) 
GAMY(I)  *  PI2INV*(FTAN*C0STHE( J+KJ)  -  FLOG*SINTHE(J+KJ)) 
GBMY(I)  «  PI2INV*(FLOG*COSTHE( J+KJ)  +  FTAN*SINTHE( J+KJ)) 
AMY(I)  «  GAMY(I)*COSALF(I)-GBMY(I)*SINALF(I) 

BMY(I)  =  GAMY( I )*SINALF( I )+GBMY( I )*COSALF( I ) 

IF  (J.EQ.NP1)  GOTO  20 
SUMAMY(I)  =  SUMAMY(I)  +  AMY(I) 

SUMBMY(I)  «  SUMBMY(I)  +  BMY(I) 

VY( I)  *  VY( I)  +  AMY(I)*QK( J+LJ) 

VX(I)  =  VX(I)  -  BMY( I)*QK( J+LJ) 

20  CONTINUE 

C 

C  ADD  THE  CONTRIBUTION  DUE  TO  THE  WAKE  ELEMENTS 
C 

VY( I)  =  VY( I)  +  SUMBMY( I )*GAMK( L)  +  SS(L)*BMY( I)* 

+  (GAMMA(L)-GAMK(L))/WAKE(L) 

VX(I)  *  VX(I)  +  SUMAMY( I)*GAMK(L)  +  SS(L)*AMY(I)* 

+  ( GAMMA( L) -GAMK( L) ) /WAKE( L) 

25  CONTINUE 

C 

C  CALCULATE  INFLUENCE  COEFFICIENT  DUE  TO  THE  WAKE 
C 

DO  35  L  =  l.NAIRFO 
KMC  =  (L-1)*MM1 
DO  30  MC  =  1.MM1 

IF  ((L.EQ.  I)  .AND.  (MC.EQ.N))  GOTO  30 
DX  =  XMID  -  XC(L,MC) 

DY  «  YMID  -  YC(L,MC) 

DIST2  =  DX*DX+DY*DY 

GCMY  '  «  -PI2INV*DX/DIST2 

GCMX  =  +PI2INV*DY/DIST2 

CMY  =  GCMY*COSALF( I)+GCMX*SINALF( I) 

CMX  =  -GCMY*SINALF( I)+GCMX*COSALF( I ) 

C 

C  ADD  CORE  VORTEX  CONTRIBUTION 
C 

VY(I)  »  VY(I)  +  CMY*CV(L,MC) 

VX(I)  «  VX(I)  +  CMX*CV(L,MC) 

30  CONTINUE 

35  CONTINUE 

C 

C  CONVECTION  VELOCITY  OF  CORE  VORTICES  AT  NEXT  TIME  STEP 
C 

CCVX(I,N)  «  VX(I) 

CCVY(I,N)  «  VY(I) 

10  CONTINUE 
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15  CONTINUE 
40  CONTINUE 
RETURN 
END 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  NEVPOSN 

TRANSFORM  THE  RESPECTIVE  FROZEN  LOCAL  COORDINATES  AND  PANEL 
ANGLES  TO  THE  GLOBAL  FRAMES  OF  REFERENCE 


XSHIFT  -  GLOBAL  X  LOCATION  OF  THE  2-ND  AIRFOIL  PIVOT  POSITION 

YSHIFT  _  GLOBAL  Y  LOCATION  OF  THE  2-ND  AIRFOIL  PIVOT  POSITION 

I TYPE  ....  FLAG  FOR  DEALING  WITH  THE  VARIOUS  COMPONENTS 

*  0  :  FOR  TRANSFORMING  AIRFOIL  COORDINATES 

*  1  :  FOR  TRANSFORMING  WAKE  ELEMENTS 

*  2  :  FOR  TRANSFORMING  MOST  RECENT  CORE  VORTEX  SHED 

*  3  :  FOR  TRANSFORMING  ALL  PREVIOUS  CORE  VORTICES 

XPRM  ....  X  GLOBAL  RELATIVE  MOVEMENT  OF  THE  PIVOTS  POSITION 

YPRM  . .  . .  Y  GLOBAL  RELATIVE  MOVEMENT  OF  THE  PIVOTS  POSITION 


SUBROUTINE  NEWPOS( ITYPE) 

COMMON  /BOD/  IFUG,NL0WER,NUPPER,N0DT0T,X(202) ,Y(202) , 

+  COSTHE( 20 1 ) , S INTHEC  20 1 ) , SS( 2 ) , NP 1 , NP2 , NP3 , NP4 , 

+  NP5 .XSHIFT, YSHIFT, NAIRF0,XI(202) ,YI(202) , 

+  C0STHL(201) ,SINTHL(201) 

COMMON  /CORV/  CV(2,200)  ,XC(2,200)  ,YC(2,200)  ,M,TD,CWX(2,200) , 
+  CWY(2,200)  ,XCI(2,200)  ,YCI(2,200) 

COMMON  /GEOM/  SINALF( 2) ,C0SALF(2) ,0MEGA(2) ,UX(2) ,UY(2) .PIVOT, 
+  XPRM, YPRM 

CHECK  FLAG  TO  ASCERTAIN  TYPE  OF  TRANSFORMATION 

IF  (  ITYPE  .EQ.  1  )  GOTO  30 

IF  (  ITYPE  .EQ.  2  )  GOTO  40 

IF  (  ITYPE  .EQ.  3  )  GOTO  45 

AIRFOIL  COORDINATES  AND  LOCAL  ANGLES  TRANSFORMATION 
DO  10  I  =  l.NPl 

X(I)  «  (XI(I)-PIV0T)*C0SALF(1)+YI(I)*SINALF(1) 

Y(I)  *  (XI(I)-PIVOT)*( -SINALF( 1))+YI(I)*C0SALF( 1) 

10  CONTINUE 

DO  20  I  ■  NP1+1,2*NP1 

X(I)  -  ( XI ( I ) -PIVOT)*COSALF( 2 )+YI ( I )*SINALF( 2 )+ 

+  XSHIFT  +  XPRM 

Y(I)  -  (XI( I) -PIVOT)*( -SINALFC 2) )+YI( I)*COSALF( 2)+ 

+  YSHIFT  +  YPRM 

20  CONTINUE 

TRANSFORM  LOCAL  ANGLES 

DO  222  L  «  1,2 
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K  *  (L-1)*(N0DT0T+1) 

DO  222  J  *  l.NODTOT 

DX  ■  X( J+l+K)  -  X( J+K) 

DY  *  Y( J+l+K)  -  Y( J+K) 

DIST  «  SQRT(DX*DX  +DY+DY) 

SINTHE( J+K)  *  DY/DIST 

222  COSTHE( J+K)  *  DX/DIST 
GOTO  50 

C 

C  WAKE  ELEMENT  TRANSFORMATION 
C 

30  I  «  NP2 

X(I)  «  (XI( I ) -PIVOT)*COSALF( 1 )+YI( I )*SINALF( 1 ) 

Y(I)  *  (XI(I)-PIVOT)*( -SINALF( l))+YI(I)*COSALF( 1) 

I  «  NP2  +  1 

X(I)  *  (XI(I)-PIV0T)*C0SALF(2)+YI(I)*SINALF(2)+ 

+  XSHIFT  +  XPRM 

Y(I)  «  (XI(I)-PIVOT)*(-SINALF(2))+YI(I)*COSALF(2)+ 

+  YSHIFT  +  YPRM 

C 

C  WAKE  ELEMENT  ANGLES  TRANSFORMATION 
C 

DO  223  L  =1,2 
K  =  (L-l)*(NODTOT+l) 

J  =  NP1 

DX  *  X(NP2+L-1)  -  X(NP1*L) 

DY  =  Y(NP2+L-1)  -  Y(NP1*L) 

DIST  *  SQRT(DX*DX  +DY+DY) 

SINTHE( J+K)  =  DY/DIST 

223  COSTHE( J+K)  *  DX/DIST 
GOTO  50 

C 

C  MOST  RECENT  CORE  VORTEX  SHED  TRANSFORMATION 
C 

40  XC( 1,M)  a  (XCI( 1,M)-PIV0T)*C0SALF( 1)+YCI( 1,M)*SINALF( 1) 

YC( 1 ,M)  a  (XCI( 1 ,M)-PIVOT)*( -SINALF( 1))+YCI( 1 ,M)*COSALF( 1) 
XC(2,M)  =  (XCI(2,M)-PIV0T)*C0SALF(2)+YCI(2,M)*SINALF(2)+ 

+  XSHIFT  +  XPRM 

YC(2,M)  a  (XCI(2,M)-PIVOT)*(-SINALF(2))+YCI(2,M)*COSALF(2)+ 
+  YSHIFT  +  YPRM 

GOTO  50 
C 

C  ALL  PREVIOUS  CORE  VORTICES  TRANFORMATION 
C 

45  DO  46  I  »  l.M-l 

XC(1,I)  =  (XCI ( 1 , I ) -PIVOT)*COSALF( 1 )+YCI ( 1 , I )*SINALF( 1 ) 
YCCl.I)  »  (XCI( 1 , I) -PIVOT)*( -SINALF( 1) )+YCI( 1 , I)*COSALF( 1) 
XC( 2,1)  «  ( XCI( 2 , I ) -PIVOT)*COSALF( 2 )+YCI( 2 , I )*SINALF( 2)+ 

+  XSHIFT  +  XPRM 

YC(2,I)  a  (XCIC2,I)-PIV0T)*(-SINALFC2))+YCI(2,I)*C0SALF(2)+ 
+  YSHIFT  +  YPRM 

46  CONTINUE 
50  CONTINUE 

RETURN 

END 
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APPENDIX  B.  EXAMPLE  INPUT  DATA  FOR  PROGRAM  USPOTF2 


5 

^rk^r^rk^rk^rk^rk^rk1^k^rk^rk^rk^rk^rk1r^rk^Hrk1rMrk^r1rk^rte^rk^rk^rk^r^rk^rk1c^^^e^Hrk^rk^rk^r^r^rle^t 

AIRFOIL  TYPE  :  8. 4%  THICKNESS  SYMMETRICAL  MISES  AIRFOIL 
AIRFOIL  COORDINATES  ARE  USER  INPUT  DATA  (IFLAG  =  1) 

NLOWER  *  25  ,  NUPPER  «  25 


01  25  25 

2  0.  0  2.  0 


1.  000000 

0.  994858 

0. 980866 

0. 958884 

0.  929536 

0. 

893455 

0. 851308 

0. 803815 

0.  751753 

0. 695948 

0.  637271 

0. 

576620 

0.514918 

0. 453098 

0.  392084 

0. 332794 

0.  276105 

0. 

222865 

0. 173861 

0. 129819 

0. 091393 

0.059146 

0.  033560 

0. 

015010 

0. 003767 

0. 000000 

0. 003767 

0. 015008 

0.  033560 

0. 

059146 

0.091393 

0. 129819 

0. 173861 

0.222865 

0.276105 

0. 

332791 

0. 392082 

0. 453095 

0.514915 

0.576617 

0. 637266 

0. 

695946 

0. 751750 

0. 803815 

0.851308 

0. 893455 

0.  929536 

0. 

958884 

0. 980866 

0.  994858 

1. 000000 

0. 000000 

-0. 000782 

-0. 002784 

-0. 005721 

-0. 009351 

-0. 

013459 

-0.017837 

-0.022285 

-0.026618 

-0.030671 

-0.  034289 

-0. 

037341 

-0.039712 

-0. 041314 

-0. 042083 

-0. 041979 

-0.  040979 

-0. 

039096 

-0.036360 

-0. 032820 

-0. 028555 

-0.023651 

-0.018220 

-0. 

012379 

-0. 006259 

0. 000000 

0. 006259 

0. 012379 

0.  018220 

0. 

023651 

0. 028555 

0. 032820 

0. 036360 

0. 039096 

0.  040979 

0. 

041979 

0. 042083 

0.041314 

0. 039712 

0.037341 

0.  034289 

0. 

030671 

0.026618 

0. 022285 

0.017837 

0. 013459 

0.  009351 

0. 

005721 

0. 002784 

0. 000782 

0. 000000 

1. 000000 

0. 994858 

0. 980866 

0.958884 

0.929536 

0. 

893455 

0.851308 

0. 803815 

0. 751753 

0. 695948 

0.637271 

0. 

576620 

0.514918 

0. 453098 

0. 392084 

0. 332794 

0.  276105 

0. 

222865 

0. 173861 

0.  129819 

0.091393 

0. 059146 

0.  033560 

0. 

015010 

0. 003767 

0. 000000 

0.003767 

0. 015008 

0.  033560 

0. 

059146 

0. 091393 

0. 129819 

0.  173861 

0. 222865 

0.  276105 

0. 

332791 

0. 392082 

0. 453095 

0.514915 

0.576617 

0.  637266 

0. 

695946 

0.  751750 

, 0. 803815 

0. 851308 

0. 893455 

0.  929536 

0. 

958884 

0. 980866 

0.  994858 

1. 000000 

0. 000000 

-0. 000782 

-0. 002784 

-0.005721 

-0.  009351 

-0. 

013459 

-0. 017837 

-0. 022285 

-0. 026618 

-0.030671 

-0.  034289 

-0. 

037341 

-0.039712 

-0.041314 

-0. 042083 

-0.041979 

-0. 040979 

-0. 

039096 

-0. 036360 

-0. 032820 

-0. 028555 

-0.023651 

-0.  018220 

-0. 

012379 

-0. 006259 

0. 000000 

0. 006259 

0.012379 

0.  018220 

0. 

023651 

0. 028555 

0. 032820 

0. 036360 

0.039096 

0.  040979 

0. 

041979 

0. 042083 

0.041314 

0.039712 

0.037341 

0.  034289 

0. 

030671 

0. 026618 

0. 022285 

0.017837 

0.013459 

0. 009351 

0. 

005721 

0. 002784 

0. 000782 

0. 000000 

0.0 

0.0 

-45. 80000 

0.00  0 

.0 

.0000 

0.0 

0.00 

0.0 

0.00  0 

.0 

.000 

5. 0000 

0.065 

0.001 

0.00  4 

.  855698  • 

-.921370 

0 

0.0 

1.  216290 


0.0 

0 
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APPENDIX  C.  EXAMPLE  OUTPUT  DATA  FROM  PROGRAM 

USPOTF2 


DATA  READ  FROM  FILE  COOE  1 
S 


AIRFOIL  TYRE  :  8.4Z  THICKNESS  SYMMETRICAL  MISES  AIRFOIL 
AIRFOIL  COORDINATES  ARE  USER  INPUT  DATA  IIFLAO  ■  1) 
NLGMER  ■  25  >  NUPPER  «  25 


I  25 

25 

2  0.0 

2.0 

1.000000 

0.994858 

0.980844 

0.958884 

0.929534 

0.893455 

0.851308 

0.803815 

0.751753 

0.495948 

0.437271 

0.574420 

0.514918 

0.453098 

0.392084 

0.332794 

0.274105 

0.222845 

0.173841 

0.129819 

0.091393 

0.059144 

0.033540 

0.015010 

0.003747 

0.000000 

0.003747 

0.015008 

0.033540 

0.059146 

0.091393 

0.129819 

0.173841 

0.222845 

0.274105 

0.332791 

0.392082 

0.453095 

0.514915 

0.574417 

0.637244 

0.695944 

0.751750 

0.803815 

0.851308 

0.893455 

0.929534 

0.958884 

0.980844 

0.994858 

1.000000 

0.000000 

-0.000782 

-0.002784 

-0.005721 

-0.009351 

-0.013459 

-0.017837 

-0.022285 

-0.024418 

-0.030471 

-0.034289 

-0.037341 

-0.039712 

-0.041314 

-0.042083 

-0.041979 

-0.040979 

-0.039094 

-0.034340 

-0.032820 

-0.028555 

-0.023451 

-0.018220 

-0.012379 

-0.004259 

0.000000 

0.004259 

0.012379 

0.018220 

0.023451 

0.028555 

0.032820 

0.034340 

0.039094 

0.040979 

0.041979 

0.042083 

0.041314 

0.039712 

0.037341 

0.034289 

0.030471 

0.024418 

0.022285 

0.017837 

0.013459 

0.009351 

0.005721 

0.002784 

0.000782 

0.000000 

I. 000000 

0.994858 

0.980844 

0.958884 

0.929534 

0.893455 

0.851308 

0.803815 

0.751753 

0.495948 

0.437271 

0.574420 

0.514918 

0.453098 

0.392084 

0.332794 

0.274105 

0.222845 

0.173841 

0.^29819 

0.091393 

0.059144 

0.033540 

0.015010 

0.003747 

0.000000 

0.003747 

0.015008 

0.033540 

0.059144 

0.091393 

0.129819 

0.173841 

0.222845 

0.274105 

0.332791 

0.392082 

0.453095 

0.514915 

0.574417 

0.437244 

0.495944 

0.751750 

0.803815 

0.851308 

0.893455 

0.929534 

0.958884 

0.980844 

0.994858 

1.000000 

0.000000 

-0.000782 

-0.002784 

-0.005721 

-0.009351 

-0.013459 

-0.017837 

-0.022285 

-0.024418 

-0.030471 

-0.034289 

-0.037341 

-0.039712 

-0.041314 

-0.042083 

-0.041979 

-0.040979 

-0.039094 

-0.034340 

-0.032820 

-0.028655 

-0.023451 

-0.018220 

-0.012379 

-0.004259 

0.000000 

0.004259 

0.012379 

0.018220 

0.023651 

0.028555 

0.032820 

0.034340 

0.039094 

0.040979 

0.041979 

0.042083 

0.041314 

0.039712 

0.037341 

0.034289 

0.030471 

0.024418 

0.022285 

0.017837 

0.013459 

0.009351 

0.005721 

0.002784 

0.000782 

0.000000 

0.000000 

0.000000-45.800003 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.000000 

0.450000 

0.025000 

0.001000 

0.000000 

4.855498 

-0.921370 

1.216290 

0 

0.000000 

0 
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AIRFOIL!  II  PERIMETER  LENGTH  ■  2.018599 
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